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I.  INTRODUCTION 

Currently,  there  is  a  need  within  the  Structures 
Directorate  for  a  simple  method  to  predict  the  thermodynamic 
properties  of  decaying  rocket  exhaust  plumes.  This  is  needed 
primarily  to  determine  the  effects  of  exhaust  plume  heating 
and  loading  on  adjacent  structures  during  missile  firings. 

In  order  to  accomplish  this,  the  method  must  describe 
accurately  the  plume  parameters  of  pressure,  temperature, and 
velocity  at  discrete  axial  and  radial  locations. 

A  simple  momentum  integral  method  was  chosen  to 
accomplish  the  above  stated  goal.  Therefore,  the  purpose  of 
this  paper  is  to  present: 

1.  the  analytical  method  chosen, 

2.  an  empirical  determination  of  the  mixing  rate 
parameter,  K 

3.  the  computer  code  that  was  developed, 

4.  and  typical  results  as  compared  to  a  similar  method 
and  measured  test  data  [1]. 

The  momentum  integral  method  used  is  based  on  the 
results  of  work  done  in  the  mid-1960's  by  Donaldson  and  Gray 
[2].  The  method  predicts  the  turbulent  mixing  and  decay  of 
axially  symmetric,  compressible,  free  jets,  where  the  mixing 
gases  are  dissimilar.  This  is  an  extension  of  previous  work 
done  by  Warren  [3]  where  the  mixing  gases  were  both  air.  In 
both  cases,  it  is  assumed  that  there  is  parallel  flow,  no 
radial  or  axial  pressure  gradients,  no  shocks,  and  an 
optimally  expanded  nozzle. 

II.  DISCUSSION 

The  integral  method  of  Donaldson  and  Gray  [2]  solves  the 
equation 


d  r  r* 

—  pu2rdr 
dxj  0 


u* 


d  fr* 

—  purdr  +  r*r* 
dx  0 


(II. l) 


This  equation  states  that  the  net  rate  of  change  in  the  axial 
direction  of  the  momentum  in  a  jet  from  the  centerline  out  to 
some  radius  r*  is  equal  to  the  change  in  momentum  due  to  mass 
addition  with  velocity  u*  at  radius  r*  plus  the  change  in 
momentum  due  to  the  shear  taking  place  at  radius  r*.  Letting 
r*  =  rs ,  where  rs  is  the  radius  at  which  the  local  velocity  is 
one  half  the  value  of  the  centerline  velocity,  and  r*  =  r5,  which 
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is  the  shear  stress  at  rs  ,  equation  (II. 1)  becomes 


d  frs  d  frs 

—  pu2rdr  =  u5  — I  p urdr  -  rsrs  (II. 2) 

dxj  0  dxj  0 

Also  letting  r*  equal  some  very  large  radius,  the  two  terms  on 
the  right  hand  side  of  equation  (II. 1)  vanish.  That  is,  u*  and 
r*  approach  zero  as  r*  approaches  infinity.  This  implies  that 
the  net  change  in  total  momentum  passing  through  a  plane  normal 
to  the  axis  is  zero  or  that  the  total  axial  momentum  does  not 
change.  This  leads  to  the  following  equation 


ao 

pu2rdr  =  const  =  /?1u12r1/2  (II. 3) 

0 


A.  Velocity  Profiles 

The  method  assumes  a  mixed  viscid-inviscid  plume  model  as 
shown  in  Figure  1.  As  can  be  seen,  the  jet  consists  of  two 
distinct  regions.  The  first  being  the  core  region,  where  there 
is  embedded  a  core  of  flow  moving  at  the  nozzle  exit  velocity. 

The  core  decreases  in  radial  extent  as  one  moves  axially  down  the 
jet  to  a  point  where  it  disappears.  Outside  the  core  turbulent 
mixing  is  occurring.  The  mixing  increases  in  radial  extent  as 
the  core  decreases  and  continues  to  increase  until  the  core 
disappears . 

The  second  region,  referred  to  as  the  developed  region, 
begins  once  the  core  has  disappeared.  In  this  region,,  the 
centerline  velocity  decreases  with  axial  distance  and  the  radial 
extent  of  the  plume  continues  to  increase. 

Two  velocity  profiles  are  needed  to  describe  the  radial 
velocity  distribution.  The  velocity  profile  for  the  core  region 
is  given  by 

u  =  ux  x  exp[->  (r2  -  r-^2 )  /  (rs  2  -  rj_2)]  r>rL  (II. 4) 

and 

u  =  Uj  r£r^ 

The  velocity  profile  for  the  developed  region  is  given  by 

u  =  uc  x  exp[-i (r/rs ) 2 ]  (II. 5) 

In  the  above  equations,  r^  refers  to  the  outer  radius  of  the 
core,  r5  is  the  radius  at  which  the  velocity  is  one  half  the 
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centerline  velocity,  ux  is  the  exit  plane  velocity,  uc  is  the 
velocity  along  the  centerline,  and  X  is  equal  to  ln2. 

B.  Turbulent  Shear  Stress 

As  was  assumed  by  Donaldson  and  Gray  [2],  the  turbulent 
shear  stress  is  given  by 

r  =  K/?bUs  (3  u/<3  r)  (II. 6) 

where  K  is  the  mixing  rate  parameter,  b  is  a  local  scale 
length,  p  is  the  local  density,  and  U?  is  an  arbitrary  local 
velocity  difference.  In  the  core  region  b  is  equal  to  r5  -  r^ 
and  Us  is  equal  to  ux/2,  which  yields 

r  =  Kp(rx  -  ri)  (u1/2)  (flu/ar)  (II. 7) 

For  the  developed  region  b  is  equal  to  r5  and  Us  is  equal  to 
uc/2,  yielding 

r  =  Kprs  (uc/2)  (3u/3r)  (II. 8) 

C.  Spreading  and  Decay 

Equations  (II. 2)  and  (II. 3)  are  the  two  basic  equations 
describing  the  flow.  Also,  expressions  have  been  developed  for 

1.  u  in  terms  of  r,  rs ,  and  r^  in  the  core  region 

[Eq. (II . 4) ] ,  and  u  in  terms  of  r,  r5 ,  and  uc  in  the  developed 

region  [Eq.(II.5)] 

2.  r  in  terms  of  K,  rs ,  and  r^  in  the  core  region 

[Eq. (II.6)],  and  r  in  terms  of  K,  r5 ,  and  uc  in  the  developed 

region  [Eq.(II.7)] 

3.  and  p  in  terms  of  u  (derived  in  Appendix  A). 

Substituting  these  expressions  into  the  two  flow  equations, 
only  three  unknowns  will  remain  in  each  region.  In  the  core 
region,  K,  rs ,  and  r^  are  unknown  and  in  the  developed  region  K, 
rs ,  and  uc  are  unknown. 

Values  of  K  have  been  determined  empirically  from 
experimental  test  data.  A  series  of  free- jet  mixing  experiments 
were  performed  by  both  Donaldson  and  Gray  [2]  and  Warren  [3],  the 
results  of  which  were  analyzed  by  use  of  the  momentum  integral 
technique  to  determine  local  turbulent  mixing  rates  and  the  local 
mixing  rate  parameter,  K.  The  results  of  these  tests  are  shown 
in  Figure  2  [2].  This  figure  gives  the  mixing  rate  parameter,  K, 
as  a  function  of  local  Mach  number,  Ms .  It  was  concluded  from 
these  tests  [2]  that  a  general  relationship  exists  at  any  given 
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axial  location  between  a  local  mixing  rate  parameter  and  the 
local  Mach  number.  It  was  also  shown  that  this  parameter  was 
independent  of  the  physical  properties  of  the  mixing  gases. 

A  determination  of  the  parameter  K  was  also  made  using  the 
computer  code  developed  and  measured  test  data  from  the  MK  56 
rocket  motor  [1].  The  curve  developed  for  K  from  this 
correlation  is  also  given  in  Figure  2.  It  differs  slightly  from 
the  results  of  Donaldson  and  Gray,  but  the  same  basic  trend  is 
apparent . 

Knowing  K  as  a  function  of  Ms  ,  only  two  unknowns  and  two 
equations  remain  for  both  the  core  and  developed  regions.  This, 
therefore, provides  a  means  to  solve  for  the  decay  and  spreading 
behavior  of  the  plume  flowfield.  The  solution  to  these  equations 
is  given  in  Appendices  B  and  C.  This  solution  allows  the 
determination  of  the  velocity  at  any  given  point  in  the  plume 
flowfield,  which  enables  the  parameters  of  total  temperature, 

Mach  number,  pitot  pressure,  and  static  temperature  to  be 
calculated  and  plotted. 

III.  COMPUTER  PROGRAM 

The  method  outlined  above  was  implemented  with  a  "C"  language 
computer  program.  The  program  included  a  slight  modification 
to  account  for  nozzle  underexpansion,  which  is  common  in  many 
Army  missiles.  A  circular  arc  curve  fit  from  the  nozzle  exit  to 
the  core  boundary  was  used,  where  the  core  boundary  is  calculated 
assuming  an  ideally  expanded  nozzle.  This  curve  fit  was  first 
developed  for  use  with  the  momentum  integral  method  by  Reardon 
and  Word  [4].  The  details  of  implementing  this  curve  fit  are 
provided  in  Appendix  D. 

The  program  requires  an  input  knowledge  of  the 

1.  ambient  temperature  and  pressure, 

2.  chamber  temperature  and  pressure, 

3.  molecular  weight  of  exhaust  gas, 

4.  specific  heat  ratio  of  exhaust  gas, 

5.  nozzle  exit  Mach  number, 

6.  nozzle  exit  radius, 

7.  and  nozzle  half  angle. 

The  program  then  calculates  local  flowfield  parameters  of  Mach 
number,  pitot  pressure,  and  total  temperature  at  discrete  axial 
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and  radial  locations.  Other  parameters  such  as  velocity, 
molecular  weight,  specific  heat,  and  the  specific  heat  ratio  are 
also  calculated  and  can  be  printed  with  slight  modifications  to 
the  program. 

The  program  is  written  in  "C"  and  a  listing  is  provided  in 
Appendix  E.  A  sample  printout  is  provided  in  Appendix  F.  Also 
included  in  this  report  is  a  plot  program  which  allows  graphical 
presentation  of  the  results.  This  plot  program  will  plot  lines 
of  constant  Mach  number,  pitot  pressure,  and  total  temperature 
within  the  plume  flowfield.  The  user  is  prompted  to  choose  one 
of  these  parameters  for  plotting.  Having  chosen  the  desired 
parameter,  the  program  automatically  chooses  ten  constant  contour 
lines  for  plotting.  Any  number  of  these  can  be  overwritten, 
giving  the  user  the  ability  to  plot  any  number  of  constant 
contour  lines  of  his  choosing  up  to  ten.  Figures  3,  4,  and  5  are 
example  plots  for  each  of  these  parameters.  A  listing  of  the 
plot  program  is  provided  in  Appendix  G. 

IV.  TYPICAL  RESULTS  AND  COMPARISONS 

A  similar  semi-empirical  method  was  developed  in  1970  by 
Piesik  [1].  In  the  presentation  of  his  method,  a  comparison  was 
made  with  measured  data  obtained  by  the  Naval  Surface  Weapons 
Center,  Dahlgren,  VA,  for  the  MK56  DTRM  [5].  The  Table  lists  the 
pertinent  rocket  motor  parameters  and  geometry  for  the  MK56  DTRM. 

Predictions  from  the  momentum  integral  method  were  compared 
to  both  the  method  of  Piesik  [1]  and  the  measured  data  for  the 
MK56  DTRM  [1].  Comparisons  were  made  for  both  the  axial  and 
radial  directions  of  the  plume  flowfield. 

Figures  6  and  7  give  a  comparison  of  the  centerline  total 
temperature  and  pitot  pressure.  As  can  be  seen,  the  momentum 
integral  method  agrees  very  well  with  the  method  of  Piesik  and 
the  measured  data. 

Radial  comparisons  of  pitot  pressure  are  presented  in  Figures 
8  thru  12.  Again,  agreement  with  the  predictions  of  Piesik  and 
the  test  data  is  quite  good  radially  from  the  centerline. 

Figures  11  and  12  actually  indicate  that  the  momentum  integral 
method  may  actually  predict  pitot  pressures  more  accurately 
than  the  method  of  Piesik  far  downstream.  As  mentioned 
previously,  in  order  to  predict  exhaust  plume  heating  or  loading, 
an  accurate  model  of  the  flowfield  plume  is  necessary.  As  is 
indicated  by  this  comparison,  the  momentum  integral  method  is  an 
effective  procedure  for  modeling  an  exhaust  plume. 
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V.  CONCLUSIONS 


From  the  MK56  DTRM  total  temperature  and  pitot  pressure 
comparison,  the  momentum  integral  method  showed  good  agreement 
with  measured  data  and  the  predictions  of  Piesik's  method  [5]. 
Therefore,  it  is  concluded  that  this  technique  is  a  valid  design 
tool  that  can  be  used  to  predict  flowfield  plume  parameters 
necessary  in  determining  heating  and  pressure  loadings  due  to 
exhaust  gas  impingement  from  rocket  motors. 
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Figure  1.  Fully  expanded  exhaust  plume  structure. 
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Figure  2 .  Empirical  determinations  of  the  mixing  parameter 
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Figure  9.  Radial  pitot  pressure  comparison  at  x  =  30  ft. 
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Figure  11.  Radial  pitot  pressure  comparison  at  x  =  50  ft, 
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Figure  12.  Radial  pitot  pressure  comparison  at  x  =  70.4  ft. 
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Specific  heat  at  constant  pressure 
Enthalpy 

Integrals  defined  in  Appendix  C 

Mixing  rate  parameter 

Molecular  Weight 

Mach  Number 

Pressure 

Radial  distance 

Nondimensional  radius,  r/r1 

Universal  gas  constant 

Temperature 

Velocity 

Nondimensional  velocity,  u/Uj 
Axial  distance 

Nondimensional  distance,  x/rx 

Specific  heat  ratio 

Nozzle  half  angle 

In  2 

Density 

Shear  stress 

Prandtl-Meyer  expansion  angle 
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Stagnation  condition 

Value  at  the  exit  of  an  underexpanded  nozzle 


Subscripts 

i 


5 

c 

d 

i 

t 

u 

* 


Exit  plane  conditions  of  fully  expanded  nozzle 

Ambient  conditions 

Value  where  u  =  luc 

Centerline  conditions 

Downstream  of  normal  shock 

At  viscid-inviscid  boundary  (core  boundary) 

Chamber  conditions 
Upstream  of  normal  shock 

Variable  value  at  some  arbitrary  radial  location,  r 
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APPENDIX  A 


DEVELOPMENT  OF  AN  EXPRESSION  FOR  DENSITY 
AS  A  FUNCTION  OF  VELOCITY 


APPENDIX  A.  Development  of  an  Expression  for  Density  as  a 

Function  of  Velocity 


In  the  following  discussion,  it  is  assumed  that  free 
mixing  is  taking  place  between  an  initially  irrotational  and 
axially  symmetric  stream  of  an  ideal  gas  given  by  the 
subscript  1  and  another  ideal  ambient  gas  given  by  subscript 
<*> .  The  gases  are  assumed  to  be  inert  chemically. 

Takinq  the  usual  Crocco  integral  (and  assuming  it  holds 
throughout  the  mixing  region) 

h  +  (u2/2 )  =  Au  +  B  (A. 1) 

and  applying  the  boundary  conditions  h  =  h^  when  u  =  0  and  h  =  hx 
when  u  =  ux  ,  the  above  equation  becomes 

h  =  hffl+(h1“-hto)(u/u1)-(u1V2)(u/u])2  (A. 2) 

where 

h: •  =  hx  +  ux  2/2  (A. 3) 

Also  taking  the  concentration  integral 

c<*  +  aau  -  b  (A-4> 

and  considering  the  mixing  of  a  rocket  exhaust  with  the  air 
surrounding  it,  there  exist  two  mass  fractions  (ca) ,  cx  and  cm , 
where  cx  is  the  local  mass  fraction  of  the  exhaust  gas  and  c^  is 
the  local  mass  fraction  of  the  ambient  or  surrounding  gas. 
Applying  the  boundary  conditions,  cx  =  0  when  u  =  0  and  cx  =  1 
when  u  =  ux  ,  equation  (A. 4)  becomes 

cx  =  u/Uj  (A. 5) 

Also  applying  the  boundary  conditions  c^  =  1  when  u  =  0  and 
cOT  =  0  when  u  =  ux  ,  equation  (A. 4)  becomes 

Coo  =  1  "  (u/ux)  (A.  6) 

Knowing  the  mass  fractions,  expressions  for  molecular  weight  and 
specific  heat  become 

m  =  1/ [  (l/mx -l/iO  (u/ux  )  +  (l/mai)  ]  (A. 7) 

Cpi  =  ( Cp x  +Cp,-o  )  (u/Uj  ) +CpOT  (A.  8) 

Assuming  h  =  CpT  with  a  constant  pressure  throughout  the  plume 
flowfield  equal  to  the  ambient  pressure,  and  using  the  perfect 


A-l 


gas  law,  density  can  be  written  as 

P  =  PmCp/(Rh)  (A. 9) 

Substituting  expressions  (A. 2),  (A. 7),  and  (A. 8)  into  (A. 9) 

yields 

Pco  [  (Cpi  +Cpco  )  (U/U1  )  ] 

p  = - 

[(l/ir^-l/m*,)  (u/Uj )  +  ( l/m^, )  ]  [hc0+(h1  °-hm)  (u/u1)-(u1 2 /2)  (u/ux)2] 

(A. 10) 

where  density  is  now  a  function  of  known  quantities  and  velocity. 
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APPENDIX  B 

DEVELOPMENT  OF  THE  EQUATIONS  FOR  THE  PLUME  FLOWFIELD 
FOR  BOTH  THE  CORE  AND  DEVELOPED  REGIONS 


APPENDIX  B.  Development  of  the  Equations  for  the  Plume 
Flowfield  for  both  the  Core  and  Developed  Regions 


Bl.  General  Equations 

The  first  step  in  solving  equations  (II. 2)  and  (II.  3)  is 
to  non-dimensionalize  these  equations  and  the  expressions  for  u, 
Eqs.(II.4  and  5),  r,  Eq.(II.8  and  9),  and  p ,  Eq.(A.lO).  Letting 
capital  letters  represent  non-dimensionalized  parameters  such 
that  R  =  r/rj ,  U  =  u/Uj ,  and  X  =  x/r1 ,  these  expressions  become 


d 

dX 


Rs  d 

(p/ P® ) U2  RdR  -  (Uc/ 2 ) 

0  dX 


fRc 


(p/p.)URdR  = 


[r5/(pcou12)  ]RX  (B.l) 


'  aa 

(p/ Poo )  U2  RdR  =  l/2(Pl/Pao)  ( B  •  2 ) 

0 

where 

U  =  XJl  x  exp[-X  (R2-Ri2)/(RS2-Ri2)  ]  (B.3) 

(t5/pw Ul2)  =  (ps/P«)[l/2Kc(R5-Ri)](3U/aR)  (B.4) 

for  the  core  region, 

U  =  Uc  x  exp[-X (R/Rs ) 2 ]  (B. 5) 

(r  5/ PajUj  2  )  =  (p  5 / Poo  )  1/2 KRS  Uc  (dU/dR)  (B.6) 

for  the  developed  region,  and 

(P/P.)  =  [ ($U+l)/(ru+l) (1+AU-BU2) ]  (B. 7) 

$  =  ( Cpi  / Cpco )  -1  r  =  (m00/m1 )  -1 


A  =  (hx  °/hc0 )  -1  B  =  l/2(u12/hco) 

B2 .  Development  of  Equations  for  the  Developed  Region 
From  Eq.(B.5)  it  can  be  shown  that 

(3U/3R)  =  -(XUc/R5)  ( B .  8 ) 

and 
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(B.9) 


RdR  =  - (Rs  2dU/2XU) 

Expressions  (B.5),  (B.6),  (B.7),  (B.8),  and  (B.9)  can  be 

substituted  into  the  two  flow  equations  (B.l)  and  (B.2). 
(Note:  For  compactness,  the  density  ratio  is  carried  as  p/pm, 
though  it  is  a  known  function  of  u  given  by  Eq.(B.7).)  These 
substitutions  and  the  necessary  change  in  the  limits  yield 


d 

ruc  i 

d 

fuc  1 

Rs2 

(P/P*) UdU 

-  (Uc/2)  — 

Rs2 

(p/p*) dU 

dX 

u5 

dX 

Us  J 

-X2KR,U  2  (p/Pao) 


(B. 10) 


and 


U, 


0 


(p/p*)  UdU  =  X/Rs2  (Pl/pJ 


(B. 11) 


Differentiating  by  parts,  Eq.(B.lO)  becomes 


CU, 


Rs  2  — 
dX 


(P  / Pea )  UdU  -  (Uc/2 )  Rs  2  — 
Uc  +  dX 


Ur 


u. 


(p/p*)  dU  + 


u. 


(P/P*)  UdU 


U. 


dR5  2 


dX 


-  (Uc/2) 


U, 


Uc 


(P/P*) dU 


dRs2 


dX 


-X2KR5Uc2  (P5/Pc0) 


(B. 12) 


The  five  indicated  integrations  in  Eqs.(B.10  and  B.ll)  can  be 
solved  closed  form  and  are  represented  by  the  following  notation 


fU, 


Ue 


(P/P*)  UdU  =  Ic  -  I, 


(B. 13 ) 


J 


U, 


Uc 


(P/P*) dU  =  Jc  -  Js 


( B . 14 ) 


d 

dX 


ruc  duc  dus 

(p/pm)UdU  =  Fc - Fe  - 


(B. 15) 


Uc 


dX 


dX 


B-2 


d  ruc  duc  dU5 

—  (P/Pco)  dU  =  Gc - G5  - 

dXJ  Us  dX  dX 


(B. 16) 


(p/ Pa ) UdU  -  Ic  -  I0 


(B. 17) 


where  I,  J,  F,  and  G  are  functions  of  U  and  are  solved  in 
Appendix  C. 

Substituting  expressions  (B.13-B.17)  in  Eqs.(B.ll)  and 
(B. 12) ,  yields 


lc  -  I0  =  \  (Pi/P co)/Rs 


(B. 18) 


— us 
2  2 


1 

duc 

-  Fs 

)  R5  2 

+ 

2 

dX 

]  dRs  2 

duc 

(Fc  ~ 

Is) 

-  = 

Jduc 

dX 

x2krsuc2 

(Ps/Pa) 

(B. 19) 


Equation  (B.18j  can  be  rearranged  to  solve  for  Rs 2  as  follows 


R  *  = 


MPx/P oo)/(Ic  -  Io) 


(B. 20) 


dR5  2 


dlc  _  R52Fc 
(Ic  -  I0)  dUc  (Ic  -  I0) 


(B. 21^ 


where 


Fc  dIc/dUc 
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Substituting  Eqs.(B.20)  and  (B.21)  into  Eq. (B.19),  provides 


i 


ruc  i  i  1 

ruc  11 

— (Gc-  — G5)-(Fc-  -Fs) 

i 

►n 

o 

— (jc-js)-(ic-is) 

*to 

to 

to 

CM 

y  X 


^2(PS/P»)UC2(IC-I0) 


dUc 

Rs  -  =  K 

dX 


(B. 22 ) 


which  is  of  the  form 


f (Uc) (dUc/dX)  =  K  (B. 23) 

Due  to  the  complex  nature  of  f( Uc)  and  the  fact  K  is  a  function 
of  Ms  and  therefore  X,  an  equation  of  this  form  must  be  solved 
numerically.  Since  it  is  a  first  order,  non-linear,  ordinary 
differential  equation^ it  was  solved  using  the  Runge-Kutta  method 
in  conjunction  with  a  Taylor  series  approximation. 


B3 .  Development  of  Equations  for  the  Core  Region 

In  the  core  region,  Eqs.(B.8)  and  (B.9)  become 

(3U/3R)  =  -ARS/(RS2  -  Rj.2)  (B.24) 

and 

RdR  =  [(R52  -  Rj_2  )  dU/2XU]  (B.25) 

As  was  done  previously,  expressions  (B.3),  (B.4),  (B.7) ,  (B.24), 

and  (B.25)  can  be  substituted  into  equations  (B.l)  and  (B.2). 
Performing  the  necessary  integration  and  using  the  integral 
notation  of  Eqs. (B. 13-B. 17) ,  two  expressions  similar  to 
Eqs.(B.18)  and  (B.19)  can  be  obtained.  (Note:  The  integration  in 
the  core  region  must  be  broken  into  two  parts.  These  are  0  <  R  < 
Rj_,  where  conditions  are  constant,  and  Rj_  £  R  <  °° . ) 


\{p,/Pao)  (1  -  Ri2) 


B-4 


*  (P1/Pa>)  (1  “  Ri2  ) 

R  2  =  -  +  Rj  2 


(B. 28 ) 


'  I0 


and 


d(R52  -  Ri2)  \{pl/pat)  dRi2 


dX 


(B. 29 ) 


I,  -  I0  dX 


Substituting  Eqs.(B.28)  and  (B.29)  into  Eq.(B.27)  and  letting 


<P  =  \(p1/pa3)/(Il  -  I0) 


(Pi/ps)  r  1 

- i - 

>  Ldr-Io) 


(J5-J1)-2(Is-I1) 


1 

J 


(B. 30) 
(B . 3 1 ) 


equation  (B.27)  becomes 

V»Ri(R5  +  Ri) 

-  dRi  =  KcdX 

<t>  +  (1  -  <p)  Ri2 


(B. 32) 


Since  il>  and  K  are  constant  in  the  core  region  Eq.(B.33)  can  be 
solved  in  closed  form,  yielding 


Ri* 


Ri2 


1  |_  [<£+  ( l-0 )  Ri 2  ]  1 


4>+  ( l-<t> )  Ri2  J 


dRi  =  K 


dX  ( B . 3  3 ) 


This  integration  will  yield  two  solutions,  depending  on  the  value 
of  0.  When  0  is  greater  than  zero  the  integration  yields 
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V>  f  [<#>+  (l-<£)  Rj_2  ]  5+R-£-2  <p 

X  =  -  i -  +  -  x 

K  L  1  -  <P  T  ( 1-0 ) 


where 

t  =  [*(wn* 


B3 .  Solving  for  Local  Flowfield  Parameters 


Having  solved  the  necessary  equations  to  determine  the 
velocity  at  any  point  within  the  flowfield,  flowfield  parameters 
of  interest  can  be  calculated  as  follows 

CP 

=  Cp^U'  +  l) 

(B. 37) 

m 

=  m*,/  (TU'+l) 

(B. 38) 

7 

=  l/(l-Rg/mCp) 

(B. 39) 

T 

=  To,  ( 1+AU-BU2  )  /  ($U+1) 

( B .  4  0 ) 

M 

=  (U)  (u1)/(7RgT/m)-^ 

( B .  4 1 ) 

P“ 

r  7-1  k 

=  P  1  +  -  M2 

( B  •  4  2 ) 

l  2  j 


B-6 


(7  +  1)Mu2 

c 

7  +  1 

_  (7  -1)  Mu2  +2j 

_27Mu2  -7  +  1. 

fp  O 


T 


7-1 

1  +  -  M2 

2 


where 


U'  =  u/uc 

C  =  7  /  ( 7  —  1 ) 
v.  =  1/  ( 7  —  1 ) 


( B .  4  3  ) 


( B . 44 ) 
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APPENDIX  C 


DERIVATION  OF  I,  J,  F,  and  G 


APPENDIX  C.  Derivation  of  I,  J,  F,  and  G 


From  Eqs.(B.13)  and  (B.17) 

I  =  J*(p/p.)UdU 

Substituting  Eq. (B.7)  in  for  p/pm  in  Eq.(C.l),  yields 

($U2+U)  dU 


I  = 


(ru+1) ( l+AU-BU2 ) 


(C.l) 


(C.2) 


This  integral  can  be  solved  in  closed  form  by  utilizing  partial 
fractions.  The  integrand  can  be  expressed  as 


*1 


£2u+£: 


(C.3) 


ru+l  (l+AU-BU2) 
where  £l#  £2,  and  £3  are  constants.  Equation  (C.2)  now  becomes 


I  - 


r  du 


ru+i 


+  a. 


UdU 


( l+AU-BU2 ) 


+ 


dU 


(l+AU-BU2 ) 


(C.4) 


Carrying  out  the  integration  yields 
£  1  £  2 

I  =  —  ln(ru+l)  -  In (l+AU-BU2)  - 

r  2B 


£2A 


B 


2£: 


(4B+A2 ) ^ 


tanh-1 


-2BU+A 


(4B+A2 ) 5 


(C.5) 


The  constants  £lf  £ 2 ,  and  £3  can  be  evaluated  by  equating  like 
terms 


3>U2+U 


£2u+e: 


(ru+i) (1+au-bu2 )  ru+i  (i+au-bu2) 


(C.6) 


such  that 


$  =  -exB  +  £2r 
i  =  etA  +  c2  +  c3r 

0  -  *1  +  ^3 


C-l 


yielding 


-($-D 


€,  = - 

( C .  7 ) 

B+rA-r2 

B+4>  (A-r) 

£2  -  - 

(C.8) 

B+rA-r2 

(*-r) 

^3  =  - - 

(C.9) 

B+rA-r2 


From  Eq. (B. 14) 

J  =  SiP/Pao)  dU  (C.10) 

Substituting  Eq. (B.7)  for  p/pal  Eq.(C.lO)  becomes 


J 


($U+l)dU 

(TU+1) (1+AU-BU2 ) 


(C. 11) 


As  above  the  integral  can  be  solved  in  closed  form  using  partial 
fractions.  Equation  (C.ll)  becomes 


'  dU 

4-  X 

UdU 

4-  C 

dU 

ru+i 

+  6  2 

+  °  3 

1+AU-BU2  J 

1+AU-BU2 

and  integrating  yields 


where 


*1 

J  =  —  ln(ru+l) 
r 


—  In ( 1+AU-BU2 )  - 
2B 


a  2A 

-  + 

L  B 


28 


1 

(4B+A2 ) i 


tanh-1 


-2BU+A 
(4B+A2 ) 1 


r  (®-D 
B+rA-r2 
B(«-n 
B+rA-r2 


(C .  14 ) 


(C. 15) 

(C. 16) 


C-2 


B+r  (a-$>) 

«s3  =  -  (C.  17 ) 

B+rA-r2 

From  Eq.(B.15)  F  can  be  defined  as  follows 
F  =  dl/dU 

which  is  identical  to  Eq.(C.3),  yielding 

F  =  -  +  -  (C .  18 ) 

ru+l  (1+AU-BU2) 

where  £2,  and  £3  have  been  previously  defined  in  Eqs. (C.7-9) . 

Also  from  Eq. (B.16),  G  can  be  defined  as 
G  =  dJ/dU 
which  leads  to 


J  =  -  +  -  ( C .  18  ) 

ru+l  (l+AU-BU2) 

where  <S  3  ,  <S  2 ,  and  <S3  have  been  defined  in  Eqs. (C. 15-17) . 
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APPENDIX  D 


MODIFICATIONS  TO  ORIGINAL  METHOD  TO  ALLOW 
FOR  NOZZLE  UNDEREXPANSION 


APPENDIX  D.  Modifications  to  Original  Method  to  Allow  for  Nozzle 

Underexpansion 


The  integral  method  outlined  in  Appendices  A,  B,  and  C 
applies  to  ideally  or  fully  expanded  nozzles.  Many  Army  missiles 
have  rocket  nozzles  that  are  slightly  underexpanded.  To  allow 
for  some  underexpansion;  the  method  was  modified  with  a  circular 
arc  curve  fit  [4].  The  details  of  this  modification  are  given 
below. 

The  radius  R  in  the  momentum  integral  method  is 
nondimensionalized  with  respect  to  rx ,  where  r1  is  the  radius  at 
the  exit  of  a  fully  expanded  nozzle.  In  most  cases;  the  rocket 
motor  to  be  analyzed  will  not  be  fully  expanded.  Given  an 
underexpanded  rocket  motor  with  a  nozzle  exit  radius  given  by  rl  ' 
and  an  exit  Mach  number  given  by  Mx ' ,  the  flow  must  be  expanded 
to  ambient  pressure  to  obtain  the  necessary  parameters  of  rx  and 
M, .  The  Mach  number  at  ambient  pressure  is  given  by 


where 


M, 


2 

r 

rpti 

V 

-  l 

7,-1 

_ 

LpcoJ 

(D.l) 


V  =  (71“1)/71 

and  the  associated  exit  radius  for  a  fully  expanded  nozzle  is 


Ml  ' 

2+(7i-1)M12 

_2+(7i-1)M1  '2. 

where 

V  =  (7l+l)/[4(7l-l)  ] 

Using  the  fully  expanded  parameters  of  rx  and  Mx ,  the  core 
region  boundary  Rj^  can  be  calculated  according  to  the  momentum 
integral  relations  developed  in  Appendix  B.  Then  using  the  real 
radius  rr '  a  circular  arc  can  be  drawn  from  the  exit  of  the 
underexpanded  nozzle  out  to  the  point  of  full  expansion  where,  r 
=  rt .  The  circular  arc  from  the  real  (underexpaned)  nozzle  and 
the  core  boundary  for  the  ideal  nozzle  will  overlap  as  is  shown 
in  Figure  D.l.  Merging  these  curves  will  create  a  smooth  crve 
describing  the  core  region,  which  expands  from  the  lip  of  the 
exit  of  the  underexpanded  nozzle  out  to  ambient  pressure  and  then 
behaves  as  per  the  momentum  integral  equations. 
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The  slope  of  the  circular  arc  at  the  exit  of  the 
underexpanded  nozzle  is  assumed  to  be  the  Prantl-Meyer  expansion 
angle  given  by 


1 

1 

u>  =  —  tan-1 

02  (Mx  2 -1) 

-  —  tan-1 

4?2  (Mt  ' 2 -1) 

_  ft. 

_  s 

where 

tan-1 

Mx  2-l 

+  tan-1 

0  = 

(7  i  ~1)/  (7  !+l) 

Mx '2-l 


+  9 


(D.  3 ) 


Knowing  the  expansion  angle,  the  axial  extent,  b,  of  the 
expansion  region  is 

b  =  tan(w) (rt '-a)  (D.4) 

and  the  radius  r  of  the  circular  arc  can  be  determined  as  a 
function  of  x,  given  by 


r  = 


(rx '-a) 2+tan2 (o) (r: '-a) 2- 


[x-ran  {i 


where 


a  =  (r1-r1 'sec(u) )/ (l-sec(w) ) 

The  local  Mach  number  can  be  calculated  from  the  area  Mach  number 
relationship 


r 

2+(7l-l)M2 

r  ' 

1  *< 

M 

_2+(7i-1)M1  '2_ 

which  provides  Mach  number  as  a  function  of  x.  This  equation  for 
Mach  number  was  solved  using  a  Newton  iterative  technique. 

Knowing  the  Mach  number,  the  velocity,  pressure,  and  temperature 
within  the  expansion  region  (x  <  b)  can  be  calculated  from  basic 
isentropic  one-dimensional  relations.  The  values  of  7,  Cp,  amd  m 
are  assumed  constant  throughout  the  core  in  this  region  and  equal 
to  their  values  at  the  nozzle  exit. 


Since,  with  the  momentum  integral  method  all  linear  values 
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are  nondimensionalized  with  respect  to  a  fully  expanded  nozzle, 
r^  and  r^  begin  at  the  ideal  nozzle  lip.  In  order  to  conform  to 
the  physical  geometry  of  an  underexpanded  nozzle,  r-j_  and  r5  must 
begin  at  the  underexpanded  nozzle  lip.  This  is  achieved  by  using 
the  relation 

Ar  =  rx  -  r 

where  r  is  given  by  Eq. (D.5) .  The  Ar  is  subtracted  from  r^  and 
r5  within  the  expansion  region  (x  <  b) ,  forcing  these  curves  to 
begin  at  the  exit  of  the  underexpanded  nozzle.  Outside  the 
expansion  region  (x  >  b)  Ar  is  equal  to  zero. 


D-3/(D-4  blank) 


APPENDIX  E 
PROGRAM  LISTING 

24  June  1988  -  01:08  PM  (Pages  E-l  thru  E-13) 
24  June  1988  -  01:13  PM  (Pages  E-14  thru  E-18) 


»»  plmrpt.c  06/24/88  01:08p 


page:  001 


/*  EXHAUST  PLUME  PROGRAM  */ 
/*  */ 
/*  Prepared  by  */ 
/*  R.U.NOURSE  */ 
/*  STRUCTURES  DIRECTORATE  V 
/*  RESEARCH,  DEVELOPMENT,  AND  ENGINEEERING  CENTER  */ 
/*  U.S.ARMY  MISSILE  COMMAND  */ 
/*  REDSTONE  ARSENAL,  AL  V 
/*  V 
/*  */ 
/*  This  program  utilizes  the  Momentum  Integral  Technique  */ 
/*  to  predict  the  thermodynamic  properties  of  decaying  */ 
/*  rocket  exhaust  plumes.  This  program  will  accurately  */ 
/*  predict  the  plume  parameters  of  pressure,  velocity,  */ 
/*  and  temperature  at  discrete  axial  and  radial  locations.  */ 

r  */ 
/*  */ 
/*  The  user  is  prompted  by  the  program  for  the  required  */ 
/*  input.  This  input  includes;  */ 
/*  */ 
/*  1.  Ambient  pressure  and  temperature  */ 
/*  2.  Chamber  pressure  and  temperature  */ 
/*  3.  Specific  heat  ratio  of  exhaust  gases  */ 
/*  4.  Molecular  wieght  of  exhaust  gases  */ 
/*  5.  Mach  number  at  the  exit  of  the  nozzle  */ 
/*  6.  Exit  nozzle  radius  */ 
/*  7.  Nozzle  half  angle  */ 
/*  */ 
/*  */ 
/*  V 
^include  <stdio.h> 

#include  <math.h> 

^include  <fcntl.h> 

double  t2,p 2;  /*  Prop  of  air  temp,  pres  */ 
double  tO, pO, garni, xmw1,xma1;  /*  Prop  of  propellant  */ 
double  rl.thel;  /*  Exit  plane  geom  */ 
double  gam2  =  (double)  1.4;  /*  Specific  heat  ratio  */ 
double  rgas  =  (double)  49720.;  /*  Real  gas  constant  */ 
double  xmu2  *  (double)  28.97;  /*  Molecular  wt  of  air  */ 
double  xk  =  (double)  0.025;  /*  Mixing  factor  */ 
double  cp1,cp2;  /*  Specific  heat  */ 
double  hi , h2;  /*  Enthalpy  */ 
double  rho0,rho1 ,rho2,rho5;  /*  Density  */ 
double  tl;  /*  Temperature  at  nozl  exit  */ 
double  al;  /*  Speed  of  sound  at  nozl  exit  */ 
double  consl  =  (double)  144.0  *  (double)  32.174; 

double  xmb;  /*  Reference  nozl  mach  number  */ 
double  rb;  /*  Reference  nozl  radius  */ 
double  omga;  /*  Expansion  angle  */ 
double  aa; 

double  xb;  /*  Axial  length  of  Region  I  */ 
double  t;  /*  Temperature  */ 
double  asound;  /*  Speed  of  Sound  */ 
double  u1;  /*  Exit  Plane  Velocity  */ 
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double  uc;  /*  Centerline  velocity  */ 
double  u5;  /*  Half  velocity  uS  =  .5  (  uc)  */ 
double  a,b,xm,p,xlam;  /*  Constants  of  integration  */ 
double  uux.uu;  /*  Non-dimensionalized  velocity  ratio  */ 

/*  U  =  u/ul  V 


double  xi ,xj ,xi 1 ,xj1 ,xi0; 
double  xi5,xj5; 
double  psi,phi,tau; 
double  xf,xg; 

double  xf5,xfc,xg5,xgc,xic,xjc; 


double  z1,z2,z3,z4,z5,z6;  /*  Parameters  for  D.E.  */ 

double  xxmax;  /*  Non-dimensionalized  max  length  of  */ 

double  xmax;  /*  Maximum  length  of  core  region  */ 

double  x;  /*  Axial  distance  */ 

double  rr,r;  /*  Radius  of  circular  arc.  Region  1  */ 

double  dr;  /*  OR  =  RB  •  R  */ 

double  xmai;  /*  Mach  nunber  Region  l  */ 

double  rri;  /*  Non-dimensionalized  radius  of  core  */ 

double  ri;  /*  Radius  of  core  region  */ 

double  rro;  /*  Non-dimensionalized  radius  of  plume  */ 

double  ro;  /*  Radius  of  plume  boundary  */ 

double  rr5;  /*  Non-dimensionalized  half  velocity  */ 

/*  radius,  RR5  =  R5  /  RB  */ 

double  r5;  /*  Half  Velocity  Radius  */ 


double  uuc.uucc; 
double  due; 
double  del.delx; 

double  cp, xmw, gam, mach , mach5 , presu , presd; 

double  xlimit,  xprt,  ent; 

int  fh; 

struct  _dat 

C 

double  dis; 
double  radii  1); 
double  tem[11]; 
double  mac  Cl  1] ; 
double  pretll]; 

>  dat; 

FILE  *prt; 

main() 

( 

int  i; 
int  byte; 
int  j; 
int  ic; 
double  temp; 
double  vC4]; 
double  uC31; 
void  boundO; 
void  calc(); 


>>>> 
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void  regionlO; 
void  viscidO; 
void  solvdeO; 
void  mix(); 
void  factorO; 
double  desubO; 
double  tstep; 

double  tstepl  3  (double)30.0; 
double  tstep2  3  (double)O.I; 
double  tstep3  =  (double)O.S; 
if  (( fh=creat( "PLUME. DAT",0666))  ==  -1) 

C 

printf ("COULD  NOT  OPEN  OATA  FILE"); 
exit(1 ); 

> 

if  ((prt  3  fopen( "PLUME. PRT»,“w"))  =3  NULL) 

C 

printf ("\n\nC0ULD  NOT  OPEN  PRINT  FILE\n\n“); 
exit(1); 

> 

pr i nt f ( "\nPROPERT I ES  OF  AMBIENT  AIR  (TEMP/PRES)  :  "); 

/*  scanf ( "XfXf ", &t2 , &p2 ) ;  */ 

t2  3  540.;  p2  3  14.7; 

pr i nt f ( " \nPROPER TIES  OF  PROPELLANT  GAS  (T0,P0,GAM1 ,XMW1 ,XMA1 )  :  "); 
scanf  (  "XfXfXfXfXf" ,  &t0 , &p0 , Sgaml , ixmw 1 , &xma 1 ) ; 
printf ("\nEXIT  PLANE  GEOMETRY  (R1.THE1)  :  ">; 
scanf < "Xf Xf » ,&r1,ithe1); 
pr i nt f ( " \nMAX I MUM  AXIAL  DISTANCE  ;  "); 
scanf ( "Xf " , &x l imi t ) ; 
printf("\nPRINT  INTERVAL  :  »); 
scanf ("Xf",&xprt); 
cnt=xprt; 

printf ("XCOUNT  3  X10.4f\n",cnt); 

/* 

PROPERTIES  OF  AM8IENT  AIR 

*/ 

dat.rad(0)=t0; 
dat.radC1)3p0; 
dat.radC2] sgaml; 
dat.rad(3]  =xmw1; 
dat.radC4]  =xma1; 
dat.rad  C5]=r1; 
dat.rad[6)=the1; 


if  ((byte  3  write(fh,(char*)&dat, (unsigned  int)sizeof(dat)))  33 


perrorC"'); 
fprintf(prt,"\n 
fprintf(prt,"\n 
fprintf (prt,"\n 
fprintf(prt,"\n 
fprintf(prt," 
fprintf (prt," 
fprintf (prt," 


PROPERTIES  OF  AMBIENT  AIR\n"); 

TEMP  s  X10.4f  PRES 

ROCKET  MOTOR  PARAMETERS'^"); 

COMBUSTION  TEMP  =  X10.4f  CHAMBER  PRES 

SPECIFIC  HEAT  RATIO  3  X10.4f  MOLECULAR  WEIGHT 

RADIUS  OF  EXIT  *  %10.4f  NOZZLE  HALF  ANGLE 

EXIT  MACH  NUMBER  =  %10.4f \n\n",xma1 ); 


•1) 


=  X10.4f\r>\n",t2,p2); 

=  X10.4f\n»,t0,p0); 

3  X10.4f\n", garni, xmwl); 
=  %10.4f\n", rl , thel ); 
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/*  Calculate  ambient  specific  heat  */ 

cp2  =  (gam2*(rgas/xmw2))/(gam2  -  1.0); 

/*  Calculate  ambient  air  enthalpy  */ 
h2  =  cp2  *  t2; 

/*  Calculate  ambient  air  density  */ 

rho2  =  (p2  *  consl  *  xmw2)  /  (rgas  *  t2); 


/* 

CALCULATE  PROPERTIES  OF  PROPELLANT  GAS 
SPECIFIC  HEAT 
V 


cpI  =  (gam1*(rgas/xmw1))/(gam1  -  1.0); 

/*  Stagnation  entalphy  */ 
hi  =  cpI  *  tO; 

/*  Density  */ 

rhoO  =  (consl  *  pO  *  xmul)  /  (rgas  *  tO); 

/*  Speed  of  sound  */ 

tl  =  tO  /  (  (double)  1.0  ♦  (double)  0.5  *  (garni  ■  (double)  1.0) 
*  (xmal  *  xmal)); 
al  =  sqrt( garni *(rgas/xmu1 )*t1 ); 
bound ( ); 

/*  Calculate  local  velocity  for  reference  nozl  */ 

t  =  tO  /  (  (double)  1.0  ♦  (double)  0.5  ‘(garni  -  (double)  1.0) 

*  (xmb*xmb)); 

asound  =  sqrt(gam1*(rgas/xmw1 )*t); 
uc  =  xmb  *  asound; 
u1*uc; 
uuc'uc/ul ; 

/*  Calculate  constants  used  for  integration  */ 

a  =  (  hi  /  h2  )  -  (double)  1.0; 
b  =  (  ul  *  ul  )  /  (  (double)  2.U  *  h2); 
xm  =  (  xmw2  /  xmwl  )  •  (double)  1.0; 
p  =  (  cpI  /  cp2  )  -  (double)  1.0; 
xlam  =  log(ldouble)  3.0); 

uux=(double)1 .0; 
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rhol  =  rho2*(p*uux+(double)1 .0)/((xnt*uux»(double)1 .0) 

*((double)1 .0+a*uux-b*uux*uux)); 

calc(O.O); 
xi  1=xi ; 
xj1=xj; 

uux=(double)0.; 
calc(0,0); 
xiO=xi ; 

u5=(double)0.5*u1  ; 
uux=u5/u1  ; 
calc(1,0); 

if  (phi  <  (double)I.O) 

( 

xxmax=(psi/xk)»((sqrt(phi  )•  (double)2.0)/(( doubled  .0-phi ) 
♦(phi/(tau*((double)1 . 0-phi )))»atan(tau/phi }); 

>  else 
{ 

xxmax=(psi/xk)*(((double)2.0-sqrt(phi ))/(phi -(double)l .0) 
+(phi/((double)2.0*tau*(phi -(double)l .0)))*log((phi -tau)/(phi-*-tau))); 

xmax=xxmax*rb; 

fprintf (prt,"  MAXIMUM  EXTENT  OF  CORE  REGION  =  %10.4f\n\n\n\n\n‘‘,xmax>; 

del  =  tstep2; 

x  =  (double)O.O; 
whi le(x  <=  x l imi t) 

< 

if  (x  <=  xb) 

< 

temp=r1-aa; 

temp*=temp; 

r*sqrt(tefrp*tan(omga)*tan(omga)*temp-  (x-  tan(omga)*(r1  -aa))* 

<x- tan(omga)*(r1 -aa) ) >+aa; 
regionl ( ); 

asound  =  sqrt(gam1*(rgas/xmw1)*t>; 
uc=xma i 'asound; 

uuc=uc/u1;  */ 

viscidO; 
dr=rb-r; 
ri=rri*rb-dr; 

rr5=sqrt(phi+( (double)l .0-phi )*rri'rri ); 

r5=rr5*rb-dr; 

tefnp=rri  -dr/rb; 

temp*=temp; 

rro=sqrt( (<(pow( rr5-dr/ rb, (doubt e)2.0) • temp)*log( (double) .02)) 
/-xlam)+temp); 
ro=rro'rb; 

>  else 

if  (x  <=  xmax) 

{ 

viscidl); 

ri=rri*rb; 

rr5=sqrt(phi+((double)1 .0-phi )*rri'rri ); 
r5=rr5'rb; 

rro=sqrt( ( ( rr5*rr5- rri'rri )*tog((double)0.02))/-xlarr»+rri*rri  ); 
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ro=rro*rb; 

mix(); 

>  else 

if  (x  >  xmax) 

( 

uc=uuc*u1 ; 

u5=(double)Q.5*uc; 

uux=u5/u1  ; 

calc(1,1); 

xi5=xi ; 

xj5=xj ; 

xf5=xf ; 

xg5=xg; 

uux=uc/u1 ; 

calc(1 , 1  ); 

xic=xi ; 

xjcsxj; 

xfc=xf; 

xgc=xg; 

rr5=sqrr(xlain*(rho1/rho2)/(xic-xi0)); 

r5=rr5*rb; 

uu=uux/uuc; 

cp=cp2* ( p*uu+ ( doub l e ) 1 . 0 ) ; 

xmw=xmw2/ <  xm*uu+ ( doub le)1.0); 

gam=(double)1.0/((double)1.0-rgas/(xmw*cp)); 

t = 1 2*  ( (  doub l e ) 1 . 0+a*uux - b*uux*uux ) / ( p*uux+ ( doub l  e )  1 . 0 )  ; 

mach=uux*u1/sqrt(gam*rgas*t/xmw); 

mach5=mach; 

factorc); 

Zl=xlam*xlam*(rho5/rho2)*(xic-xi0); 

z2=xgc- (double). 5*xg5; 

z3=xfc- (double) .5*xf5; 

z4=xjc-xj5; 

z5=xic-xi5; 

z6=xic-xi0; 

ic=0; 

v[0]=(double)0.0; 

v(1]  =uuc; 

v(2]s(double)0.0; 

v[3]=(double)0.0; 

sol vde(v,u, del,. 05,1, &ic); 

duc=vt2]  ; 

delx=del/rb; 

uucv=duc*delx; 

rro=rr5*sqrt( log( (double) .02)/-xlam); 

ro=rro*rb; 

mix(); 

) 

if  (x  <  xmax) 

( 

if  (fabs(x-xmax)  <  0.1) 

x=xmax; 

else 
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x+=xmax/(double)60.0; 

>else 

( 

if  (x  <  xmax*(double)1 .2) 

C 

x+=tstep2; 

det=tstep2; 

}else 

C 

x+=tstep3; 

del=tstep3; 

> 

> 

> 

close  (fh); 
fclose(prt); 

> 

double  desub(v) 
double  v[] ; 

extern  double  z1,z2,z3,z4,z5,z6,rr5,xlc,xfc; 
return  (xk/rr5)*((vm*vm*z1)/(z6*((vM]/(double)2.0)*z2-z3) 
-xfc*((vCl]/(double)2.0)*z4-z5))); 

> 

void  bound ( ) 

{ 

double  b0,b1 ,b2,b3(b4(b5,b6,b7,b8; 
double  rbO.rbl; 
double  bet; 

/* 

Calculate  reference  nozl  mach  number 

*/ 

b0=gam1+<double)1 .0; 
b1=gam1  •  (double)  1.0; 

xmb=sqrt(  (  (double)  2.0/bl)  *  (pow(  p0/p2,  bl/gaml)  -  (double)  1.0)); 

b2=b0/(  (double)  4.0*b1); 

b3=xmb*xmb  •  (double)  1.0; 

b4=xma1*xma1  -  (double)  1.0; 

bet=sqrt(b1/b0); 

rb0=  (double)  2.0+bl  *  xmb*xmb; 

rb1=  (double)  2.0+bl  *  xmal+xmal; 

rb=sqrt(xma1/xmb)*pow(  rbO/rbl,  b2)  *  rl; 

b5=  (double)  1.0/bet; 

b8=bet*bet; 

b6=sqrt(b8*b3); 

b7=sqrt(b8*b4); 

C*nga=b5*(atan(b6)-atan(b7))-atan(b3)+atan(b4)+the1  ; 

aa=(rb-r1»(  (double)  1 .0/cos(omga) ))/(  (double)  1.0-  (double)  1 .0/cos(omga) ); 
xb=tan(omga)*(r1 -aa); 

> 

void  calc(opt1 ,opt2) 

int  optl;  /*  opt  =  0  initial  conditions  */ 
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int  opt2; 

double  c1,c2; 
double  e1,e2,e3; 
double  g1,g2,?3; 

double  alpha, beta, gamma, delta, epsn.zeta; 

double  v1,v2,v3; 

cl  =  b  ♦  a  *  xm  -  xm  *  xm; 

c2  =  sqrt  (  a  *  a  +  (double)  4.0  *  b  ); 

el  =  (  b  *  (  p  -xm  ))  /  cl; 

e2  =  xm*((b  +  p*(a-  xm))  /cl  ); 

e3  =  (  xm  *  (  p  -  xm  ))  /  cl; 

gl  =  (  b  *  xm  *  (xm  -  p>)  /  cl; 

g2  =  -gl; 

g3  =  xm  *  ((  xm  *  (  a  -  p  )  +  b)  /  cl  ); 

alpha  =  -el  /  (  b  *  xm  ); 

beta  *  -e2  /  (  (double)  2.0  *  b  *  xm  ); 

gamma  =  (  e3  /  xm  +  (  a*  e2)  /  ((double)  2.0  *  b  *  xm  )) 

*  ((double)l .0/c2); 

delta  *  -g1/(b*xm); 

epsn  *  -g2/((double)2.0*b*xm); 

zeta  =  (g3/xm*(a*g2)/((double)2.0*b*xm))*((double)1 .0/c2); 

vl  =  xm*uux+(double)1 .0; 

v2  =  (double)I.O  +  a*uux  -  b*uux*uux; 

v3  =  (c2+(double)2.0*b*uux-a)/(c2-(double)2.0*b*uux+a); 

xi  =  alpha*log(v1 )+beta*log(v2)+gamna*log(v3); 
xj  »  delta*log(v1)+epsn*log(v2)+zeta*log(v3); 

if  (optl  *=  0) 
return; 
if  (opt2  ==  0) 

C 

xi5»xi; 

xj5=xj; 

rho5=rho2*(p*uux+(double)1 .0)/((xm*uux+(double)1 .0) 

*( (doubl e) 1 . 0*a*uux ■ b*uux*uux )); 
psi=((rho1/rho5)/xlam)*(((double)1 .0/(xi 1-xi0))*((xj5-xj1) 
-(double)2.0*(xi5-xi 1))- (doubt e)1 .0); 
phi=xlam*(rho1/rho2)/(xi1-xi0); 

if  (phi  >  (double)l.O)  tau=sqrt(phi*(phi-(double)1.0)); 
else  tau=sqrt(phi*((double)1 .0-phi )); 

return; 

> 

xf=( -el )/(b*v1 )+(e2*uux+e3)/(xm*v2); 
xg=( -gl )/(b*v1 )+(g2*uux+g3)/(xm*v2) ; 
rho5=rho2*(p*uux+(double)1 ,0)/((xm*uux+(double)1 .0) 

* ( ( doub  t  e ) 1 . 0+ a  *  uux - b*uux  *uux ) ) ; 


> 

void  regionK) 

< 
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double  xo; 
double  dfo; 
double  xn; 

double  tol  =  (double)l .e-6; 
double  tolx; 
double  test; 
double  d1,d2,d3,d4; 
int  j; 

/* 

This  subroutine  utilizes  Newtons  method  to  solve  for  the  mach 
number  in  Region  l  at  a  specific  location  x  less  than  b. 

V 

xo=xma1 ; 

d1=xma1/((r/r1 )*(r/r1 )); 
d3=gam1 - (double)l .0; 
d2=(double)2.0+d3*xma1*xma1; 
d4=(gam1+(double)1 ,0)/((double)4.0*d3); 
for  ( j=0; j<20; j++) 

< 

dfo=(double)2.*d4*d1*pow(((double)2.i-d3*xo*xo)/d2, 

(double)2.*d4- (double) 1 . )*((Cdouble)2.*d3*xo)/d2)-(double)1 . 
xn=xo-(d1*pow(((double)2.*d3*xo*xo)/d2, (double)2.*d4)-xo)/dfo; 
test=fabs(xn-xo); 
tolx=tol*fabs(xn); 
xmai=xn; 

if  (test  <  tolx) 

< 

return; 

> 

xo=xn; 

> 

> 

void  viscidO 

double  xo; 
double  dfo; 
double  xx ; 
double  xn; 
double  test; 
double  tolx; 

double  tol  =  (double)l .e-4; 
double  w1,w2,w3,w4,w5,w6,w7,w8, w9; 
int  j; 

/* 

This  subroutine  calculates  the  core  boundary  using  Newton's  method. 

*/ 

xo=(double)1.0; 
xx=x/rb; 
if  (phi  >  1.0) 

< 

for( j=0; j  <20 ; j++) 
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< 

ul =( (xk*xx)/ps i )*(phi- (double) 1 .0); 
w2=sqrt(phi -(phi- (double) 1 .0)*xo*xo); 
w3=phi/((double)2.*tau); 
w4=(phi+tau*xo)*(phi-tau); 
w5=(phi-tau*xo)*(phi+tau); 

w6=(double) .5*pow(phi ■ (phi -(double) 1 .0)*xo*xo, (double) -0.5); 
w7=(double)-2.0*(phi -(double)l .0)*xo; 
w8=tau*(phi -tau); 
w9=(-tau)*(phi+tau); 

dfo=(-w6)*w7- (double)1.0+w3*(w8/w4-w9/w5); 

xn=xo- ( (doubt e)2. -w2-xo+w3*( log(w4)- log(w5))-w1 )/dfo; 

test=fabs(xn-xo); 

tolx=tol*fabs(xn); 

rri=xn; 

if  (test  <=  tolx) 
return; 

xo=xn; 

> 

>else 

( 

for(J=0; j<20; j++) 

{ 

w1=((xk*xx)/psi )*( (double) 1.0- phi ); 
w2=sqrt(phi+((double)1 .0-phi )*xo*xo); 
w3=phi/tau; 
w4=atan(tau/phi ); 
w5=atan((xo*tau)/phi ); 

w6=(double)0.5*pow(phi+((double)1.0-phi)*xo*xo,(double)-0.5); 

w7=(double)2.0*((doubte)1 .0-phi )*xo; 

w8=(double)1 .0/((double)1 .0+(xo*tau)/phi ); 

dfo=w6*w7+(double)1 .0-w8; 

xn=xo-(w2+xo-(double)2.0+w3*(w4-w5)-wl )/dfo; 

test=f abs(xn- xo) ; 

tolx=tol*fabs(xn); 

rri=xn; 

if  (test  <=  tolx) 
return; 
xo=xn; 

> 

> 


void  mix( ) 

( 

double  temp; 
int  i; 
int  byte; 
double  mac2; 
dat.dis=x; 
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if  (<x  >  xb}  £S  (x  <  xmax)) 

< 

temp=(ro-ri )/(double)10.; 
dat. rad  101 =ri  -  temp; 

)else 

C 

temp=ro/(double)10.; 
dat.rad[0)=(double)0.0  •  temp; 

> 

for  ( i=0; i <10; i++) 

< 

if  (i  >  0)  dat.radCi]  =  dat.radCi-1]  +  temp; 
else  dat.rad[03  ♦=  temp; 
if  ((x  >  xb)  &&  (x  <  xmax)) 

C 

rr=dat.rad[i]/rb; 

uux=uuc*exp( -xlam*(rr*rr- rri'rri )/(rr5*rr5-rri*rri )); 

>else 

{ 

rr=dat.rad£i]/rb; 

uux=uuc*exp( -xtam*((rr*rr)/(rr5*rr5))); 

> 

uu=uux/uuc; 

cp=cp2*(p*uu+(double)1 .0); 

xmw=xmw2/(xm*uu+(double)1 .0); 

gam=(double)1 .0/((double)1 .0-rgas/(xmw*cp)); 

t  a  1 2*(  ( doub  t  e )  1 . 0+a*uux  -  b*uux*uux )  /  ( p*uux+  ( doub  le)1 .0);' 

dat.mac[i]=uux*uI/sqrtCgam*rgas*t/xmw); 

mac2=dat.mac(i]*dat.macti3 ; 

dat. tem[i]=t*((double)1 .0*( (gam- (double) 1 ,0)/(double)2.0)*mac2); 
presu=p2*pow( (double) 1 .0+( (gam- (doubt e)1 . 0)/( doubt e)2.0)* 
mac2,gam/(gam-(double)1 .0)); 
if  (dat.macCi)  <  (double)I.O) 

( 

dat . pre [ i ] =presu; 

}else 

( 

dat .pre[i] =presu*pow(( (gam+(double)l .0)*mac2)/ 

((gam- (double)l .0)*mac2+(doubte)2.0) , gam/ (gam- (double) 1 .0)) 

*pow((gam+(double)1.0)/((double)2.0*gam*mac2-gam«-(double)1 .0), 

(double)l .0/(gam- (double)l .0)); 

> 

dat.preti] -=p2; 

) 

/*  if  ((x  <  xmax  )  ||  (x  >  20.))  */ 

if  ((byte  =  write(fh,(char*)&dat, (unsigned  int)sizeof(dat)))  ==  -1) 
perrorC1"); 
if  (dat.dis  >=  xprt) 

C 

fprintf(prt, "AXIAL  0IST  '•); 
fprintf (prt,"%10.4f\n",dat.dis); 
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> 

void 

C 


>>>> 


fprintfCprt, "RADIUS 
for  Ci=0;i<10;i*+) 

fprintf Cprt, "XI 0.4f",dat. radii] ); 
fprintfCprt, "W); 
fprintfCprt, "MACH  NO.  '•); 

for  Ci=0;i<10;i++) 

fprintfCprt, "XI 0.4f",dat.mac[i] ); 
fprintfCprt, "Xn"); 
fprintfCprt, "TOTAL  TEMP  “>; 

for  Ci=0;i<10;i++) 

fprintf Cprt,"X10.4f",dat.temti] ); 
fprintfCprt, "\n"); 
fprintfCprt, "TOTAL  PRES 

for  C i=0; i <10; i++) 

fprintfCprt, "XI 0.4f",dat.pre[i] ); 
f  pr i n  t  f  C  pr t , "\n\n" ) ; 

xprt+  a  cnt; 

> 


factorO 


int  m; 
int  n; 

double  deplS]; 
double  indl5]; 

dep [0] = C doub t e )0 . 042 ; 
deptl] aCdouble)0. 042; 
dep[2]acdouble)0.03; 
depC3]=Cdouble)0.025; 
dep[4] =Cdouble)0.025; 

ind[0] =Cdouble)0.0; 
ind£1]=Cdouble)0.3; 
indC2]=tdouble)0.73; 
ind[3]aCdouble)1.0; 
indC4]=Cdouble) 100.0; 


if  CCmachS  >  indtO])  &&  Cmach5  <=  indll])) 
C 


m*1 ;  n=0; 

> 

if  CCmachS  >  indll])  &&  CmachS  <=  ind[2])) 
C 


m=2;  n=1; 

> 

if  CCmachS  >  indt2])  &&  Cmach5  <=  ind[3])> 
C 


m=3;  n=2; 

> 

if  CCmachS  >  indt3])  &&  Cmach5  <=  indl4])) 
C 


m=4;  n=3; 
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> 

xk=dep£n]  *  ((depCnti  -  depCn])/(  ind[m]  •  ind[n]  ))*(mach5  -  indtn]); 
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/*  This  program  solves  a  differential  equation  using  */ 
/*  using  a  Runge-Kutta  solution.  */ 
^include  <stdio.h> 

(Kinclude  <math.h> 


#def i ne  TRUE  1; 
#def ine  FALSE  0; 
double  desubO; 


void  sotvde(v,w,del ,delmin,n, ic) 
double  v£] ; 

double  wd;  /*  3  component  working  array  */ 

double  del;  /*  Integration  step  size  */ 

double  delmin;  /*  Min  value  to  use  when  decreasing  step  size  */ 

int  n;  /*  Order  of  equation  to  solve  */ 

int  *ic;  /*  Used  for  internal  control  */ 


< 

int  DBG  =  0; 
int  npl; 
int  np2; 
int  nh; 
int  loopl; 
int  loop2; 
int  i , j ; 
douMe  tstop; 
double  h; 
double  h2; 
double  vs  £22); 
double  c  [4]  £201; 
double  f; 
double  tmpl; 
double  tmp2; 


/*  INITIALIZE  */ 

npl  =  n+1; 
np2  =  n+2; 

/*  COMPUTE  STOP  FOR  THIS  INTEGRATION  STEP  */ 

tstop  *  v£0]+del; 

if  (DBG) 

printf ("NP1  =  Xd  NP2  =  %d  TSTOP  =  %lf\n“,np1 ,np2, tstop); 

/*  CHECK  FOR  FIRST  TIME  INTO  ROUTINE  */ 

if  (DBG) 

printf  ("DEL  =  %lf  W[1]  =  r.lf  IC  =  XdVV'.del.wEI]  ,*ic); 

if  (*ic  <=  0)  /*  FIRST  TIME  IN  */ 

( 

/*  RESTART  IF  INPUT  DEL  HAS  CHANGED  */ 
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if  (del  !  =  will) 

( 

*ic  =  0; 
h  =  del; 
nh  =  0; 
w  Cl  1  =  del; 

>  else 
C 

h  =  h  £0] ; 

nh  =  (int)wtl); 

/*  TEST  FOR  DOUBLING  H.  INCREASE  TEST  AS  NUMBER  OF  HALVINGS 

INCREASES.  */ 

if  (nh  !=  0) 

C 

if  (*ic  >=  (10  *  (nh  +  1))} 

< 

*ic  ♦=  1; 
h2  =  2.0  *  h; 

if  (modf ((tstop  +  0.01  *  h  -  v[0]),  &h2)  < 
C 

h  =  h2; 
nh  -=  1; 


> 

/*  SET  LOOP  COUNTERS  TO  TRUE  */ 

loopl  =  TRUE; 
loop2  *  TRUE; 

/*  INCREASE  COUNTER  */ 

while  (loopl) 

{ 

if  (DBG) 

{ 

printf  ("DEL  =  %lf  U Cl]  *  Xlf  IC  *  ltd  H  =  Xlf\n»,del  ,wCD  ,*ic,h); 

printf  ("V[0]  =  Xlf  VC1]  =  Xlf  V C2]  =  XI f  V[3]  =  Xlf\n»,vt01  ,v[1]  ,vC2]  ,vC33  ); 

) 

*ic  +=  1;  /*  label  60  */ 

/*  INTEGRATE  USING  R-K 

SAVE  V  TABLE  VALUES  IN  VS  ARRAY  */ 

for  ( i=0; i <np2; i++)  /*  label  100  */ 

vs[i]  =  v C i ] ; 

/*  FIRST  PASS  THRU  R-K  */ 

while  (l oop2 ) 

( 
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for  (j=0; j<n; j++)  /*  label  110  */ 

{ 

if  (j  ==  (n- 1 ) ) 

C 

f  =  desub(v); 

>  else 
{ 

f  *  vtj+2] ; 

> 

c CO]  Cj]  =  f  *  h; 

if  (DBG) 

printf ("1  F  =  XI f  CCOHXd]  =  Xlf\n“,f.  j.cCO]  Cj] ); 

> 


/* 


SECOND  AND  THIRD  PASS  THRU  R-K  */ 

vCOJ  =  vsCO]  +  0.5  *  h; 

for  (i*1 ; i<3; i++)  /*  do  129  */ 

{ 

for  (j=1;j<np1; j++)  /*  do  122 

< 

v  £  j  3  =  vsCj]  ♦  0.5  *  c  C  i  - 1  ]  Cj  - 13 ; 


> 

for  <  j=0;  j<n;j++)  /*  do  127 

< 

if  (j  -=  In- 1)) 

C 

f  =  desub(v); 

}  else 
t 

f  =  vtj+23 ; 

3 


cCU  Cjl  =  f  *  h; 

if  (DBG) 

printf  ("2  F  =  Xlf  CCXdlCXd]  =  Xlf\n",f ,  i ,  j  .cCUCj]  ); 

} 

if  (DBG) 

printf  ("V CO]  =  Xlf  VC1]  =  Xlf  VC2]  =  Xlf  VC3]  =  Xlf\n‘\vC0]  ,vC1]  ,vC2]  ,vC3] ); 
) 


if  (DBG) 

C 

printf (MDEl  =  Xlf  MCI]  =  Xlf  IC  =  Xd  H  =  Xlf\n‘\del , wCU  ,*ic,h); 
printfC‘0.5  *  H  =  Xlf  DELMIN  =  Xlf\n»,0.5*h,delmin); 

> 

if  ((0.5  *  h)  >=  delmin) 

C 

/*  TEST  FOR  HALVING  H  */ 


*7 


V 


tmp2  =  0.0; 

for  (j=0;j<n;j++)  /*  do  153  */ 

( 

tmpl  =  c CO]  Cj]  -  cCU  Cj]; 

if  (tmpl  !=  0.0)  tmp2  =  (cCUCj]  -  c C2]  C j] )  /  tmpl 
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if  (D8G) 

printf ("TMP2  *  Xlf\n",tmp2); 


if  (fabs(tmp2)  >  0.025) 


C 


>  else 

> 


h  =  0.5  *  h; 
nh  +*  1 ; 

*ic  =  1; 

for  ( j=0; j<np2; j++) 
vtj] 

break; 


loop2  =  FALSE; 


/*  label  155 


vsCj]; 


> 

if  (DBG) 

printf ("1/2  H  DEL  =  Xlf  W Cl  1  =  Xlf  1C  »  %d  H  =  Xlf\n",del,wCD ,*ic,h); 
)  else 
t 


loop2  = 

FALSE; 

> 

>  /*  whi le  ( loop2) 

*/ 

FOURTH  PASS  THRU  R-K 

*/ 

vCO]  =  vs  CO]  +  h; 

/*  label 

160 

*/ 

for  ( j=1; j<np1; j++) 

C 

vCj]  =  vsCj]  + 

cC2]  Cj-D; 

/*  label 

161 

*/ 

> 

if  (DBG) 

printf ("4th  pass  VCO]  =  Xlf  VC1]  =  Xlf  V[2]  =  XI f  V[3]  =  Xlf\n",vC0]  ,vC1)  ,v[2] ,  vC3] ); 
for  ( j=0; j<n; j++) 

C 

if  (j  ==  (n-1)) 

( 

f  *  desub(v); 

>  else 

f  =  vlj+2]; 

} 

c[3)  tj)  =  f  *  h; 

> 


/*  UPDATE  V  TABLE  */ 

for  ( j=0; j<n; j++)  /*  label  180  */ 

( 

vCj+D  =  vslj+1]  +  (cCO)Cj)  +2.0  *  (cCD  tj]+cC2Hj] )  +  cC3]Cj])/6.0; 

> 

if  (DBG) 

printf  ("UPDATE  VCO]  =  Xlf  VCD  =  Xlf  VC2]  =  Xlf  VC33  =  Xlf\n",vC0]  ,vC1] ,  vC21  ,vC3]  ); 

/*  TEST  FOR  END  OF  INTEGRATION  STEP  */ 
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/*  label  200  */ 

if  < C v tOl  +  0.1*h)  >  tstop) 

< 

loopl  =  FALSE; 

>  else 
C 

if  (*ic  >=  (10  *  (nh  +  1))) 

( 

*ic  +=  1; 
h2  =  2.0  *  h; 

if  (modf ((tstop  ♦  0.01  *  h  •  v[0]>,  &h2)  <  (0.011  *  h)) 

C 

h  =  h2; 
nh  - =  1  ; 

) 

> 

loop2  =  TRUE; 

) 

) 

/»  ENO  OF  REQUIRED  INTEGRATION  •/ 


v[0]  =  tstop;  /*  label  220  */ 

v[np2  -  1]  =  desub(v); 
w [01  =  h; 
w  Cl  1  =  del; 
w(2]  =  (double)nh; 


> 


return; 
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EXAMPLE  PROGRAM  PRINTOUT 


1 
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PROPERTIES  OF  AMBIENT  AIR 

TEMP  =  540.0000  PRES  =  14.7000 


ROCKET  MOTOR  PARAMETERS 

COMBUSTION  TEMP  =  4500.0000  CHAMBER  PRES  =  1004.0000 

SPECIFIC  HEAT  RATIO  =  1.2300  MOLECULAR  WEIGHT  =  28.0900 

RADIUS  OF  EXIT  =  j.5000  NOZZLE  HALF  ANGLE  =  0.2620 

EXIT  MACH  NUMBER  =  2.7800 

MAXIMUM  EXTENT  OF  CORE  REGION  =  149.7483 


AXIAL  D 1ST 

RADIUS 

102.328 

4.666 

5.883 

7.100 

8.317 

9.535 

10.752 

11.969 

13.186 

14.403 

15.620 

MACH  NO. 

3.234 

2.535 

1.965 

1.507 

1.142 

0.851 

0.618 

0.434 

0.290 

0.183 

TOTAL  TEMP 

4499.997 

4031.254 

3497.186 

2933.653 

2382.740 

1884.600 

1468.362 

1146.530 

915.468 

760.668 

TOTAL  PRES 

172.074 

104.157 

60.616 

33.542 

17.301 

8.521 

4.198 

1.991 

0.875 

0.346 

AXIAL  01  ST 

RADIUS 

200.250 

0.000 

2.669 

5.339 

8.  oca 

10.677 

13.347 

16.016 

18.686 

21.355 

24.025 

MACH  NO. 

2.347 

2.249 

1.990 

1.648 

1.293 

0.969 

0.692 

0.467 

0.292 

0.166 

TOTAL  TEMP 

3724.767  3647.135 

3417.333 

3049.904 

2581.146 

2071.737 

1592.597 

1199.945 

917.533 

737.549 

TOTAL  PRES 

87.231 

79.709 

61.362 

40.494 

23.085 

11.429 

5.316 

2.306 

0.884 

0.284 

AXIAL  DIST 

RADIUS 

300.250 

0.000 

3.681 

7.362 

11.044 

14.725 

18.407 

22.088 

25.769 

29.451 

33.132 

MACH  NO. 

1.617 

1.560 

1.406 

1.191 

0.953 

0.723 

0.518 

0.346 

0.212 

0.118 

TOTAL  TEMP 

2906.058 

2844.216 

2662.110 

2374.272 

2013.681 

1630.873 

1279.786 

998.713 

800.299 

675.496 

TOTAL  PRES 

37.906 

34.973 

27.516 

18.505 

10.775 

5.752 

2.823 

1.236 

0.463 

0.143 

AXIAL  DIST 

RAO  I  US 

400.250 

0.000 

4.745 

9.491 

14.237 

18.983 

23.729 

28.475 

33.221 

37.967 

42.713 

MACH  NO. 

1.231 

1.191 

1.082 

0.925 

0.746 

0.567 

0.404 

0.266 

0.161 

0.088 

TOTAL  TEMP 

2377.585 

2325.544 

2173.988 

1939.045 

1651.627 

1353.817 

1086.588 

876.391 

729.902 

638.535 

TOTAL  PRES 

19.624 

18.136 

14.334 

9.828 

6.041 

3.360 

1.672 

0.724 

0.264 

0.079 

AXIAL  DIST 

RADIUS 

500.250 

0.000 

5.930 

11.860 

17.790 

23.720 

29.651 

35.581 

41.511 

47.441 

53.371 

MACH  NO. 

0.980 

0.950 

0.866 

0.743 

0.600 

0.455 

0.321 

0.209 

0.124 

0.067 

TOTAL  TEMP 

1993.479 

1949.517 

1822.758 

1629.594 

1398.033 

1162.822 

955.334 

794.256 

683.023 

614.050 

TOTAL  PRES 

11.050 

10.262 

8.275 

5.878 

3.722 

2.101 

1.043 

0.444 

0.158 

0.046 

AXIAL  0 1ST 

RADIUS 

O00.250 

0.000 

7.239 

14.479 

21.719 

28.959 

36.199 

43.439 

50.679 

57.919 

65.159 

MACH  NO. 

0.804 

0.780 

0.712 

0.611 

0.492 

0.371 

0.260 

0.167 

0.098 

0.052 

TOTAL  TEMP 

1708.905 

1671.755 

1565.462 

1405.626 

1216.998 

1028.281 

863.918 

737.541 

650.842 

597.302 

TOTAL  PRES 

6.882 

6.431 

5.265 

3.806 

2.439 

1.380 

0.678 

0.283 

0.098 

0.023 
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AXIAL  DIST 

700.250 

RADIUS 

0.000 

8.663 

17.326 

25.989 

MACH  NO. 

0.676 

0.655 

0.598 

0.513 

TOTAL  TEMP 

1496.904 

1465.362 

1375.637 

1242.050 

TOTAL  PRES 

4.634 

4.344 

3.583 

2.609 

AXIAL  DIST 

800.250 

RADIUS 

0.000 

10.190 

20.380 

30.570 

MACH  NO. 

0.578 

0.561 

0.511 

0.438 

TOTAL  TEMP 

1337.040 

1310.050 

1  133.601 

1120.626 

TOTAL  PRES 

3.295 

3.093 

2.559 

1.868 

AXIAL  DIST 

900.250 

RAO  I US 

0.000 

11.808 

23.617 

35.426 

MACH  NO. 

0.503 

0.488 

0.444 

0.379 

TOTAL  TEMP 

1214.804 

1191 .486 

1125.653 

1028.910 

TOTAL  PRES 

2.440 

2.291 

1.898 

1.385 
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34.653 

43.316 

51.979 

60.642 

69.306 

77.969 

0.412 

0.309 

0.214 

0.137 

0.079 

0.042 

1086.240 

932.107 

799.128 

697.596 

628.271 

585.584 

1.677 

0.945 

0.459 

0.188 

0.064 

0.0183 

40.760 

50.951 

61.141 

71.331 

81.521 

91.711 

0.350 

0.261 

0.180 

0.114 

0.065 

0.034 

990.009 

861.886 

752.122 

668.749 

612.020 

577.164 

1.199 

0.670 

0.321 

0.130 

0.043 

0.012 

47.234 

59.043 

70.852 

82.660 

94.469 

106.278 

0.302 

0.224 

0.153 

0.096 

0.055 

0.029 

917.802 

809.513 

717.230 

647.410 

600.025 

570.956 

0.885 

0.491 

0.233 

0.093 

0.031 

0.008 
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#include  <gl.h> 

#include  <device.h> 

#include  <stdio.h> 

#include  <fcntl.h> 

^include  <math.h> 

#include  <time.h> 

int  xmin  =  1 00; 
int  ymin  =  100; 
int  xmax  =  900; 
int  ymax  *  600; 
int  choice; 
long  stat; 
short  val; 
double  maxv; 
double  minv; 
double  maxh; 
double  minh; 
double  tmp; 
double  incv.inch; 

int  DEBUG  =  FALSE;  /*  Change  to  TRUE  for  DEBUG  printout  */ 

struct  _dat 
C 

double  dis; 
double  radtll]; 
double  temlll]; 
double  macCll]; 
double  pretll]; 

>; 

struct  _dat  dat; 
struct  _dat  in; 

FILE  *plt; 

FILE  *prt; 

mainO 

C 

int  fh; 
int  ch; 
int  i , j ; 
int  inc; 
int  rbyte; 

double  maxval  =  (double)-10000000.0; 
double  minval  =  (double)  1 0000000.0; 
double  xmax  =  (double)- 10000000.0; 
double  xmin  =  (double)10000000.0; 
double  ymax  =  (double) -10000000.0; 
double  ymin  =  (double)10000000.0; 
double  nval; 
double  incr; 
int  nun; 
char  tmp(803  ; 
char  stC3); 
double  x,r,z[40]  ; 


>>»  ptumeplt.c  06/24/88  02:05p 
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void  plpltO; 
ginitO; 
cursoffO; 
qdevice(M0USE1 ); 
qdevice(KEY8D); 

/*  V 

co l or (SLUE); 

clearO; 

color(YELLOW); 

if  ((fh=open("PLUME. DAT", (0_CREAT|0_RDUR), (0666)))  ==  -1) 

printf ("COULD  NOT  OPEN  DATA  FILE"); 
exit(1); 

} 

if  ((pit  =  fopen( "PLUME . PLT", "wt" ) )  ==  NULL) 

C 

pr i ntf ( "\n\nC0ULD  NOT  OPEN  PRINT  F ILE\n\n") ; 
exit(1); 

> 

if  (DEBUG) 

{ 

prt  =  f open( "PLUMEPLT .PRT" , "wt" ) ; 

> 

choice  =  0; 

while  ((choice  !  =  ' T * )  &&  (choice  !  =  ’P1)  &&  (choice  !  =  ’M1)) 

C 

cmov2i(100,740); 

charstr("Enter  Item  to  plot  <  T)emp,  P)ress,  M)ach  >  "); 

stat  »  qread(Sval); 

tmpCO]  =  val ; 

tmpCI]  a  1  \0 * ; 

charstr(tmp); 

choice  =  toupper(val ); 

> 

cmov2i (100,720); 
charstrC’Choice  =  "); 
charstr(tmp); 

rbyte  =  read(fh,(char*)&in,sizeof(in)); 

while  ((rbyte  =  read(fh,(char*)&dat,sizeof(dat)))  ==  sizeof(dat)) 
C 

swi tch(choice) 

case  1 T1  :  for  ( i =0; i < 1 0 ; i++) 

if  (dat.tem(i)  >  maxval)  maxval  =  dat.tem(i); 
if  (dat.tem(i)  <  minval)  minval  =  dat.tem(i); 
if  (dat.rad(i)  >  ymax)  ymax  =  dat.  radii); 
if  (dat.rad(i)  <  ymin)  ymin  =  dat.rad(i); 

> 

break; 

case  'P1  :  for  ( i =0; i <10; i++) 

( 

if  (dat.preli)  >  maxval)  maxval  =  dat.preli]; 
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if  (dat.preti]  <  minval)  minval  =  dat.preti]; 
if  (dat.radti]  >  ymax)  ymax  =  dat.radti]; 
if  (dat.radti]  <  ymin)  ymin  =  dat.radti]; 

> 

break; 

case  ‘M1  :  for  ( i =0; i <10; i++) 

C 

if  (dat.macti]  >  maxval)  maxval  =  dat.macti]; 
if  (dat.macti]  <  minval)  minval  =  dat.macti]; 
if  (dat.radti)  >  ymax)  ymax  =  dat.radti]; 
if  (dat.radti]  <  ymin)  ymin  =  dat.radti]; 

) 

break; 

) 

if  (dat.dis  >  xmax)  xmax  =  dat.dis; 
if  (dat.dis  <  xmin)  xmin  =  dat.dis; 

> 

swi tch(choice) 

< 

case  1 T '  ;  maxval  =  (double)((long)maxval); 

minval  =  (double)((long)minval); 
break; 

case  'P‘  :  maxval  =  ((double)(long)(maxval*10))  /10.0; 

minval  =  ((double)(tong)(minval*10))  /l 0.0; 
break ; 

case  'H*  :  maxval  =  ((double)( long)(maxvat*100))  /100.0; 

minval  =  ((double)(long)(minval*100))  /100.0; 
break; 

> 

cmov2i(100,700); 

charstrC'Min  Value  =  "); 

spr i nt  f ( tmp, "XI 2 . 4f " ( mi nva l ) ; 

charstr(tmp); 

cmov2i(100,680); 

charstrC'Max  Value  =  "); 

sprintf(tmp,"X12.4f"(maxval); 

charstr(tmp); 

inc  =  10; 

incr  =  (maxval  -  minval); 
nun  =  0; 

for  (i=0;i<10;i++) 

zti]  =  maxval  •  incr*log10((double)(i+1)); 
while  (nun  >  -1) 

< 

color(BLUE); 

rectfi(100, 100, 1000, 680); 

color(YELLOU); 

for  (i=0;i<10;i++) 

< 

cmov2 i ( 1 00 , 660 • 20* i ) ; 
charstrC'Inc  "); 
sprintf(tmp,"Xd", i ); 
charstr(tmp); 
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sprintf  (tmpf,,X12.2f",z[i]  ); 
charstr(tmp); 

> 

cmov2 i(100, 400 ) ; 

charstrC'Enter  #  to  change  or  ' -'  to  exit  ->  _ 

cmov2i (450,400); 

stat  =  qread(Sval); 

num  =  vat  -  48; 

spr i nt  f ( tmp, "%d" , num) ; 

tmpCi]  =  1  \0 1 ; 

charstr(tmp); 

if  (vat  >  47) 

C 

cmov2i(100,380); 
charstrC'Enter  VALUE  for  #''); 
charstr(tmp); 

charstrC _ _ 

cmov2i (280,380); 
i  =  0; 
do 
( 

stat  =  qread(Sval); 
if  (vat  !=  13) 

C 

tmpCU  *  vat; 

sprintf (st, "%c", vat); 
st  (1]  =  1  \0 1 ; 
charstr(st); 

> 

}  while  (vat  !=  13); 
tmpCi]  =  1  \0  •  ; 
zCnun]  =  atof  (tmp); 
if  (DEBUG) 

< 

fprintf (prt,"tmp  =  Xs  nun  =  Xd  zCnum)  =  %10.2f\n", tmp, nun, z [nun] ) ; 

} 

> 

color(BLUE); 

clearO; 

minvat  =  z[9]; 
maxvat  =  z[0]; 

f pr i nt f ( pi t , "XI 0 . 2f XI 0 . 2f \n" , ymi n, ymax) ; 
fprintf(pl t,"X10.2fX10.2f\n",xmin,xmax); 
fprintf  (pit,  "X10.2fX10.2fX10. 2f\n»,- 1.0,  -1.0,  -1.0); 
for  ( j=0;j<inc; j++) 

( 

lseek(fh, ( long)0,0); 
vat  =  FALSE; 
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while  (Crbyte  =  read(fh,(char*)&dat,sizeof(dat))>  ==  sizeof(dat)) 

( 

swi tch(choice)  ^ 

{ 

case  ‘T1  :  for  ( i=0; i <9; i»+> 

{ 

if  (((zCj]  >=  dat.temCi] )  && 

(zCj]  <=  dat.temti+1] ))  || 

((zCj]  <=  dat.temCi])  && 

(zCj]  >=  dat.temCi+1] ))) 

C 

r  =  dat. radii]  ♦ 

<(zlj]  -  dat.temCi]  )/(dat.tem(i+1]  • 
dat.temCi]])*  (dat.radCi+1]  -  dat.radCi]); 
fprintf  Cptt("X10.2ft]0.2f%10.2f\n,,(  r,dat.dis,zCj]  ); 
val  =  TRUE; 
break; 

> 

> 

break; 

case  1 P ‘  :  for  ( i =0 ; i <9 ; i ++ ) 
f 

if  C((zCj]  >=  dat.preCi] )  && 

(zCj]  <=  dat.preli+1] ))  || 

CCzCj]  <=  dat.preCi] )  && 

(zCj]  >=  dat.preCi+1]  ))> 

< 

r  =  dat.radCi]  + 

CCzCj]  ‘  dat.preCi] )/(dat.preli+1]  • 
dat.preTi]  ))*Cdat.radCi+1]  •  dat.radCi]); 
fprintf  (pit,  "%10.2f%10. 2f  5(10. 2f\n“, 

r.dat.dis.zCj] ); 

val  =  TRUE; 
break; 

> 

> 

break; 

case  • M ■  :  for  ( i=0; i<9; i++) 

( 

if  C((zCj]  >=  dat.macCi])  && 

(ZCj]  <=  dat.mac(i+1] ))  || 

CCzCj]  <=  dat.macCi])  && 

(zCj]  >=  dat.macCi*1] ))) 

( 

r  =  dat.radCi]  + 

( ( z C j]  -  dat.macCi]  )/(dat.macCi  +  1]  • 
dat.macCi]  ))*(dat.radCi+1]  -  dat.radCi]); 
fprintf  (plt,"X10.2fX10.2f5C10.2f  \n", 
r.dat.dis.zCj] ); 

val  =  TRUE; 
break; 

> 

) 

break; 
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> 

}  /*  whi le  */ 

if  (vaD  fprintf(plt,-X10.2f*10.2fX10.2f\n-, 

>  /*  for  (j  */ 

fclose(plt); 

close(fh); 

plpltO; 

stat  =  qread(Sval); 
fclose(plt); 
if  (DEBUG) 

( 

close(prt); 

> 

color(BLACK) ; 

clear(); 

cursonO; 

gexitO; 

exit(O); 


void  plpltO 
C 

int  i , j; 

int  inc; 

int  rbyte; 

int  val; 

int  i tmp; 

double  rtmp; 

double  incr; 

double  x,r,z; 

double  asave , bsa ve ,  zsave; 

int  a,b,c; 

char  str(20]; 

float  locC3); 

float  wid.hgt.rot; 

struct  tm  *tm; 

long  elock; 

time(Sclock); 

tm  *  localtime(Scloclc); 

if  ((pit  *  f open( "PLUME . PLT" , "r " ) )  ==  NUU) 

< 

pr i ntf ( "\n\nC0ULD  NOT  OPEN  PRINT  FILE\n\n“); 
exit(1); 

> 

rbyte  =  fscanf(plt,"%f%f\n",Sminv,5(naxv); 
if  (OEBUG) 

i  fprintf (prt,"Rbyte  =  %d  Min  V  =  %10.2f  Max  V  =  %10.2f\n», rbyte, minv.maxv); 
} 

rbyte  =  fscanf(plt,“%W\n^minh,&maxh); 
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if  (DEBUG) 

fprintf(prt,"Rbyte  =  Xd  Min  H  =  X10.2f  Max  K  =  X10.2f\n‘,,rbyte,minh,maxh); 

> 

/*  axis  */ 

rtmp  =  maxv  -  minv; 

z  =  (double)((int)log10((double)rtmp)); 
incv  =  pow(10.0,z); 

if  (DEBUG) 

C 

fprintf(prt,"rtmp  %12.4f  incv  X12.4f\n", rtmp, incv); 

> 

rtmp  =  maxh  -  minh; 

z  =  (double)((int)log10((double)rtmp)); 
inch  *  pow(10.0,z); 

if  (DEBUG) 

( 

fprintf(prt,"rtmp  X12.4f  inch  X12.4f\n", rtmp, inch); 

> 

if  ((minv  -  incv)  <  0.0)  minv  =  0.0; 

else  minv  =  (double)  ((int)(minv  -  incv)); 
maxv  =  ((double)  ((int)((maxv  +  incv)/incv)))  *  incv; 

if  ((minh  •  inch)  «  0.0)  minh  =  0.0; 

else  minh  =  (double)  ((int)(minh  -  inch)); 
maxh  =  ((double)  ((int)((maxh  +  inch)/inch>))  *  inch; 
if  (DEBUG) 

( 

fprintf (prt/'minv  X10.2f  maxv  X10.2f  minh  X10.2f  maxh  X10.2f\n", minv, maxv, minh, maxh); 

> 

color(BLACK); 

clearO; 

color(WHITE); 

/*  Draw  and  label  x-axis  with  grid  */ 

move2i (xmin.ymin); 
draw2 i ( xmax , ymi n ) ; 
rtmp  =  minh; 

while  (rtmp  <=  maxh* inch/2.0) 

( 

a  =  xmin  ♦  (int)(((double)(xmax  •  xmin)/(maxh  ■  minh))  *  (rtmp  -  minh)); 
move2i(a,ymin); 
draw2i (a.ymax); 
if  (DEBUG) 

< 

fprintf(prt,"a  =  Xd  ymin  Xd  ymax  Xd  rtmp  X10.2f\n",a,ymin,ymax,rtmp); 

> 
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sprintf (str,"%6.2f",rtmp); 
if  (DEBUG) 

C 

fprintf(prt,"str  %s\n",str); 

> 

wid  =  strwidth(str); 

hgt  =  getheightC'K '); 

loctO]  =  (float)(a  •  wid/2.0); 

loc(1]  =  (float)(ymin  •  hgt  ■  5.0); 

loc[2]  =  0.0; 

gennote(wid,hgt, loc,0.Q,str); 
rtmp  +=  inch; 


/*  Draw  and  label  y-axis  with  grid  */ 

move2i (xmin.ymin); 
drau2i (xmin,ymax); 
rtmp  =  minv; 
whi le  (rtmp  <=  maxv) 

C 

b  =  ymin  -  ( int)(((double)(ymin  •  ymax)/(maxv  -  minv))  *  (rtmp  -  minv)) 

move2i(xmin,b); 

draw2i (xmax.b); 

sprintf (str,"X6.2f", rtmp); 

wid  =  strwidth(str); 

hgt  =  getheightC'M"); 

toe (01  *  (float)(xmin  •  wid-10.0); 

loc[1]  =  (floatXb); 

loc(2]  =  0.0; 

gennote(wid,hgt, loc,0.0,str); 
rtmp  +a  inev; 

> 

/*  Annotate  axis  */ 

switch(choice) 

C 


case 

•T1  : 

strcpy(str, "TOTAL  TEMPERATURE"); 
break; 

case 

•  P  •  : 

:  strcpy(str, "PITOT  PRESSURE  (PSIG)"); 
break; 

case 

'M'  : 

:  strcpy(str,"MACH  NO."); 
break; 

> 

wid  =  strwidth(str); 
hgt  =  getheightC'M"); 

loc(0]  =  (f loat)((xmax  +  xmin)/2  -  ( int)(wid/2. )); 
locCI]  =  730.0; 
loc (2)  a  0.0; 

gennote(wid,hgt, loc,0.0,str); 

/*  Plot  Title  */ 
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strcpy(str, "PLUME  PROGRAM"); 

wid  =  strHidth(str); 

hgt  =  getheightC'M"); 

locCO]  =  (f loat)(xmax+xmin)/2  -  wid/2.0; 

loc Cl]  =  750; 

l  oc  C2]  =  0.0; 

gennote(wid,hgt, loc,0.0,str); 

strcpy(str, "AXIAL  DISTANCE  (IN.)"); 

wid  =  strwidth(str); 

hgt  =  getheightO'M"); 

toctO]  =  (f loat)(xmax+xmin)/2  -  wid/2.0; 

toe  Cl]  =  50; 

loc[21  =  0.0; 

gennote(wid,hgt,loc,0.0,str); 

strcpyCstr, "RADIAL  DISTANCE  (IN.)"); 

wid  =  strwidth(str); 

hgt  =  getheightC'M") ; 

locCI]  =  (f loat)(ymax+ymin)/2  •  wid/2.0; 

locCO]  =  20; 

loc  [2]  =  0.0; 

gennote(wid,hgt, loc,90.0,str); 

/*  Input  Data  */ 

hgt  =  gethe i gh  t ( "M" )  * .  8; 
strcpyCstr, "ROCKET  MOTOR  PARAMETERS"); 

Hid  =  strwidth(str); 

toe  CO]  =  (f  loat)(xfnax+xmin)/2  •  wid/2.0; 

loc  Cl]  =  690; 

locC2]  *  0.0; 

gennote(wid,hgt, loc,Q.0,str); 

strcpyCstr, "TO  Combustion  Temp.,  Deg  R  :"); 

wid  =  strwie'th(str); 

locCO]  »  (float)xmin  +  50.0; 

loc  Cl]  =  675; 

loc  [2]  =■  0.0; 

gennote(wid,hgt, loc,0.0,str); 
sprintf(str,"%10.2f", in. rad(0] ); 
wid  ■  strwidth(str); 
gennote(wid,hgt, loc,0.0,str); 

strcpy(str,"P0  Chamber  Pressure,  psi  :"); 

wid  =  strwidth(str); 

locCO]  =  (float)xmin  +  50.0; 

loc Cl]  =  660; 

loc  [2]  =  0.0; 

gennote(wid,hgt, loc,0.0,str); 
sprintf(str,"X10.2f", in.radtl] ); 
wid  =  strwidth(str); 
gennote(uid,hgt, loc,0.0,str); 

strcpyCstr, "GAM1  Specific  Heat  Ratio  :"); 

uid  =  strwidth(str); 
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locCO]  =  (float)xmin  +  50.0; 
loctl]  =  645; 
loc  [2]  =  0.0; 

germote(wid,hgt, loc,0.0,str); 
sprintf(str,"X10.2f", in. rad £2] ); 

Hid  =  strwidth(str); 
gennote(wid,hgt, loc,0.0,str); 

strcpy(str,"XMU1  Molecular  Weight 

wid  strwidth(str); 

locCOj  =  (float)xmin  +  50.0; 

loc  til  =  630; 

loc C2]  =  0.0; 

gennote(wid,hgt, loc,0.0,str); 
sprintf (str,"%10.2f“, in. rad £31 ); 
wid  =  strwidth(str); 
gennote(wid,hgt, loc,0.0,str); 

strcpy(str,"XMAl  Exit  Mach  Number 
wid  =  strwidth(str); 

locCO]  a  (floatKxmax  +  xmin)/2.Q  +  60.0; 
toc£1]  =  675; 
loc £2]  =  0.0; 

gennote(wid,hgt, loc,0.Q,str); 
sprintf (str,"X8.2f“, in.rad [4] ); 
wid  =  strwidth(str); 
gennote(wid,hgt, loc,0.0,str); 

strcpy(str,"R1  Exit  Radius,  in.  :•<); 
wid  *  strwidth(str); 

loc £0]  =  (floatKxmax  *  xmin)/2.0  +  60.0; 
loc  £1]  =  660; 
loc  £2]  =  0.0; 

gennote(wid,hgt,loc,0.0,str); 
sprintf <str,"X8.2f" , in. rad £5] ) ; 
wid  a  strwidth(str); 
gennote(wid,hgt, loc,0.0,str); 

strcpy(str, "THE1  Nozzle  Half  Angle 
wid  =  strwidth(str); 

loctO]  =  (floatKxmax  ♦  xmin)/2.0  +  60.0; 
loctl]  =  645; 
loc  £2]  =  0.0; 

gennote(wid,hgt, loc,0.0,str); 
sprintf(str,MX8.2f", in. rad £6] ); 
wid  =  strwidth(str); 
gennote(wid,hgt, loc,Q.0,str); 

strcpy(str,"0ate  of  Run 
wid  =  strwidth(str); 

locCO]  =  (floatKxmax  +  xmin)/2.0  +  60.0; 
loctl]  =  630; 
loc £2]  =  0.0; 

gennote(wid,hgt, loc,0.0,str); 
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sprintf (str, 11  %02d/%02d/£02d" , tm->tm_mon+1 , tm->tm_mday, tm->tm_year); 

wid  =  strwidth(str); 
gennote(wid,hgt, l oc, 0.0, str); 

move2i(xmin+10, 710); 
draw2i(xmax- 10, 710); 
draw2i (xmax- 10,620 ); 
draw2i (xmin+10,620); 
draw2 i (xmi n+ 1 0 , 71 0 ) ; 

zsave  =  -9999.99; 

white  ((rbyte  =  fscanf(plt,“%f%f%f\n",&r,£x,&z))  !=  EOF) 

C 

a  =  xmin  ♦  ( int )( ( (double)(xmax  -  xmin)/(maxh  -  minh))  *  (x  -  minh)); 
b  =  ymin  -  (int)(((doubte)(ymin  •  ymax)/(maxv  •  minv))  *  (r  -  minv)); 

if  (x  <  0.0) 

< 

if  (zsave  •=  -9999.99) 

{ 

move2 i ( asave, bsa ve ) ; 
sprintf (str, "%6.2f“, zsave); 
wid  =  strwidth(str)*.8; 
hgt  =  getheight("a")*.8; 
toe  Cl]  ~  bsave; 
locCO]  =  asave; 
toe  [2]  =  0.0; 

gennote(wid,hgt, toe, 0.0, str); 

> 

rbyte  =  fscanf(plt,“7.fy.fy,f\n“,&r,&x,&z); 
a  =  xmin  +  ( int)(((double)(xmax  -  xmin)/(maxh  -  minh))  * 

(x  -  minh)); 

b  =  ymin  -  ( int )( ( (doubleKymin  •  ymax)/(maxv  -  minv))  * 

(r  -  minv)); 


movs2i (a,b); 

>  else 
{ 

draw2i (a,b); 
zsave  =  z; 
asave  =  a+10; 
bsave  =  b; 

> 

>  /*  whi te  */ 

return; 
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/*  OOC 

file:  gennote.c 

functions:  gennoteO,  CharStrO,  StrUidthO,  FontO,  GetHeightO 

These  functions  provide  a  set  of  primitive  vector  text  routines  which  are 
used  throughout  wsmgr  in  place  of  different  bitmapped  fonts.  Vector  text 
has  several  advantages: 

1)  It  is  available  (only  three  bitmapped  fonts  on  the  Si l  Gr  now) 

2)  It  can  be  rotated,  translated,  etc 

3)  When  ported,  the  text  size  can  be  scaled  to  the  new  resolution. 
The  routines  can  be  slow.  In  fact  a  profile  of  wsmgr  shows  that  gennote 

is  one  of  the  major  hogs  in  the  program. 

In  practice,  you  should  never  have  to  call  gennote.  CharStrO  and  the  other 
guys  provide  a  more  natural  method  of  handling  text.  They  were  designed  to 
mirror  the  bitmapped  function  calls.  See  charstrO,  strwidthO,  getheightO 
as  examples.  NOTE  that  you  should  carefully  check  the  type  and  number  of 
arguments  passed  to  these  functions.  These  functions  do  differ  from  their 
bitmap  counterparts  in  that  respect. 

StrUidthO  and  GetHeight  return  an  estimate  of  the  the  number  of  pixels 
required  to  print  a  string  in  the  current  font. 

gennoteO  was  written  by  Gregg  Geis  for  other  purposes.  The  code  was  ported 
to  the  Silicon  Graphics  and  modified  to  get  around  the  limitations  we  found 
with  bitmapped  fonts. 

END  DOC  */ 


^include  <math.h> 

^include  <gl.h> 

int  letterwidth=10,  letterheight=12,  strokewidth=1; 


/*  Load 

starting  index 

within 

vector 

array  for  each  character.  */ 

short 

chstart[128]  = 

<  977, 

987, 

1015,  1027,  1039,  1047,  1053, 

1065, 

1081, 

1097, 

1113,  1131,  1137,  1143,  1169, 

1183, 

1199, 

1215, 

1231,  1247,  1277,  1313,  1327, 

1355, 

1391, 

1411, 

1435,  1465,  1491,  1519,  1523, 

1539, 

0, 

559, 

583,  601,  623,  655,  683, 

707, 

719, 

731, 

743,  759,  769,  781,  785, 

795, 

367, 

387, 

397,  419,  445,  459,  477, 

497, 

507, 

539, 

799,  823,  849,  855,  865, 

871, 

901 , 

1, 

11,  35,  53,  67,  81, 

93, 

115, 

127, 

139,  155,  167,  175,  185, 

193, 

211, 

225, 

249,  267,  293,  301,  315, 

321, 

331, 

341, 

351,  943,  951,  955,  963, 

969, 

971 , 

1561, 

1587,  1607,  1625,  1645,  1665, 

1683, 

1709, 

1723, 

1745,  1769,  1781,  1791,  1815, 

1829, 

1847 

1867, 

1887,  1899,  1925,  1941,  1955, 

1961 , 

1975, 

1985, 

2005,  0,  0,  0,  0, 

0  ); 

/*  Load 

terminating  index  for 

character  vectors.  */ 

short 

chstopC128I  = 

<  986, 

1014, 

1026,  1038,  1046,  1054,  1064, 

1080, 

1096, 

1112, 

1130,  1136,  1142,  1168,  1182, 

1198, 

1214, 

1230, 

1246,  1276,  1312,  1326,  1354, 

1390, 
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1410, 

1434, 

1464, 

1490, 

1518,  1522, 

1538, 

1560, 

0, 

580, 

600, 

622, 

654, 

682, 

706, 

718, 

730, 

742, 

758, 

768, 

780, 

784, 

794, 

798, 

386, 

396, 

418, 

444, 

458, 

476, 

496, 

506, 

538, 

558, 

820, 

846, 

854, 

864, 

870, 

900, 

942, 

10, 

34, 

52, 

66, 

80. 

92, 

114, 

126, 

138, 

154, 

166, 

174, 

184, 

192, 

210, 

224, 

248, 

266, 

292, 

300, 

314, 

320, 

330, 

340, 

350, 

366, 

950, 

954, 

962, 

968, 

972, 

976, 

1586, 

1606, 

1624, 

1644,  1664, 

1682, 

1708 

1722, 

1744, 

1768, 

1780, 

1790,  1814, 

1828, 

1846 

1866, 

1886, 

1898, 

1924, 

1940, 

1954, 

1960, 

1974 

1984, 

2004, 

2012, 

0, 

0, 

0, 

0, 

o  >; 

/*  load  character  vector  data. 

The  character  set  data  defines  relative  vectors  within  a  14  by  22  grid  listed 
as  a  sequence  of  (x,y)  pairs.  The  normal  b ase_line  is  at  a  y-value  of  4. 

Desenders  for  the  lower  case  alphabet  can  then  extend  down  to  a  y-value  of  0. 

The  characters  are  printed  in  fixed  width  16  by  22  spaces,  so  interletter 
spacing  does  not  need  to  be  included  in  the  character  vector  data. 

An  (x,y)  with  the  values  50,50  wilt  force  a  pen-up  move  to  the  succeeding 
location,  thereby  allowing  characters  with  disjoint  components;  such  as, 
semi-colons,  lower  case  'i1,  etc.,  to  be  generated.  */ 
short  charset  [4000]  =  { 

/*  A  1..10  */  0,3,6,22,12,3,10,10,2,10, 

/*  B  11. .34  */  0,4,0,22,8,22,12,19,12,14,9,13,0,13,9,13,12,12,12,7,8,4,0,4, 

/*  C  35. .52  V  12,6,10,4,2,4,0,6,0.20,2,22,10,22,12,20,12,18, 

/*  D  53.-66  */  0,4,0,22,8,22,12,17,12,9,10,4,0,4, 

/*  E  67. .80  */  12,4,0,4,0,14,8,14,0,14,0,22,12,22, 

/*  F  81. .92  V  0,4,0,14,8,14,0,14,0,22,12,22, 

/*  G  93. .114  */  6,12,12,12,12,6,10,4,2,4,0,6,0,20,2,22,10.22,12,20,12,18, 

/*  H  115. .126  V  0,4,0,22,0,12,12,12,12,22,12,4, 

/*  I  127.. 138  V  4,4,8,4,6,4,6,22,4,22,3,22, 

/*  J  139. .154  */  0,8,0,6,2,4,6,4,8,6,8,22,6,22,10,22, 

/*  K  155. .166  */  0,4,0,22,0,12,10,22,2,14,12,4, 

/*  l  167. .174  */  0,22,0,4,12,4,12,6, 

/*  M  175. .184  */  0,3,0,22,6,12,12,22,12,3, 

/*  N  185, .192  */  0,3,0,22,12,3,12,22, 

/*  0  193. .210  */  0,7,0,19,3,22,9,22,12,19,12,7,9,4,3,4,0,7, 

/*  P  211. .224  •/  0,3,0,22,10,22,12,20,12,14,10,12,0,12, 

/*  0  225. .248  */  0, 7,0, 19  3,22,9,22,12,19,12,7,9,4,12,1,6,8,9,4,3,4,0,7, 

/*  R  249. .266  */  0,3,0,22,10,22,12,20,12,14,10,12,0,12,6,12,12,3, 

/*  S  267. .292  */  0,6,2,4,10,4,12,6,12,10,10,12,4,14,0,16,0,20,2,22,10,22,12,20,12,18, 
/*  T  293. .300  */  6,3,6,22,0,22,12,22, 

/*  U  301. .314  */  0,22,0,6,2,4,10,4,12,6,12,4,12,22, 

/*  V  315. .320  V  0,22,6,4,12,22, 

/*  U  321. .330  */  0,22,3,4,6,17,9,4,12,22, 

/*  X  331. .340  */  0,4,12,22,6,13,12,4,0,22, 

/*  Y  341. .350  V  0,22,6,14,12,22,6,14,6,4, 

/*  Z  351. .366  V  0,22,12,22,6,13,4,13,8,13,6,13,0,4,12,4, 

/*  0  367..386  */  0,6,0,20,2,22,8,22,10,20,0,6,2,4,8,4,10,6,10,20, 

/*  1  387. .396  */  4,4,8,4,6,4,6,22,4,18, 

/*  2  397. .418  */  0,18,0,20,2,22,10,22,12,20,12,16,10,14,2,12,0,10,0,4,12,4, 

/*  3  419. .444  */  0,20,2,22,10,22,12,20,12,16,8,14,2,14,8,14,12,10,12,6,10,4,2,4,0,6, 
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/*  4  445. .458  */  10,22,0,12,0,10,12,10,10,10,10,20,10,4, 

/*  5  459. .476  •/  0,8,2,4,10,4,12,6,12,12,10,14,0,14,0,22,12,22, 

/*  6  477.. 496  */  4,22,0,12,0,7,2,4,8,4,10,6,12,8,10,12,5,13,0,11, 

/*  7  497.. 506  */  0,20,1,22,12,22,3,4,3,4, 

/*  8  507.. 538  */  2,14,0,16,0,20,2,22,10,22,12,20,12,14,10,12,2,14,0,12,0,6,2,4,10,4,12,6,12,10,10,12, 
/*  9  539.. 558  */  12,14,10,12,2,12,0,14,0,20,2,22,10,22,12,20,12,13,8,4, 

/*  no  vector  necessary  for  space  character  */ 

/*  !  559. .580  */  5,4,7,4,7,6,5,6,5,4,50,50,6,8,8,20,6,22,4,20,6,8,0,0, 

/*  »  581. .600  •/  2,22,3,20,4,22,2,22,50,50,6,22,7,20,8,22,6,22, 

/*  #  601. .622  */  0,12,11,12,50,50,1,18,12,18,50,50,2,10,4,20,50,50,8,10,10,20, 

/»  S  623. .654  */  0,6,2,4,10,4,12,6,12,10,10,12,2,12,0,14,0,18,2,20,6,20,6,22,6,2,6,20,10,20,12,18, 

/*  X  655. .682  */  0,22,4,22,4,18,0,18,0,22,50,50,8,4,12,4,12,8,8,8,8,4,50,50,0.4,12,22, 

/*  &  683. .706  */  12,6,2,18,2,20,4,22,6,22,8,20,0,12,0,6,2,4,8,4,10,6,12,10, 

/*  1  707. .718  */  6,18,8,20,8,22,6,22,6,20,8,20, 

/*  f  719. .730  */  6,22,4,20,2,16,2,10,4,6,6,4, 

/*  )  731. .742  */  6,22,8,20,10,16,10,10,8,6,6,4, 

/*  *  743. .758  */  0,6,12,18,50,50,0,18,12,6,50,50,6,6,6,18, 

/*  +  759. .768  */  2,13,10,13,50,50,6,7,6,19, 

/*  ,  769. .780  V  6, 2, 8, 4, 8, 6, 6, 6, 6, 4, 8, 4, 

/*  -  781.. 784  */  3,13,9,13, 

/*  .  785. .794  «/  8, 4, 8, 6, 6, 6, 6, 4, 8, 4, 

/*  /  795. .798  */  1,4,11,22, 

/*  :  799. .820  */  5,6,7,6,7,8,5,8,5,6,50,50,5,16,7,16,7,18,5,18,5,16,0,0, 

/*;  821. .846  */  5,2,7,4,7,6,5,6,5,4,7,4,50,50,5,16,7,16,7,18,5,18,5,16,0,0, 

/«  <  847. .854  */  10,17,0,11,10,5, 

/*  =  855.-864  */  2,14,10,14,50,50,2,8,10,8, 

/*  >  865. .870  */  0,17,10,11,0,5, 

/*  ?  871. .900  */  2,18,2,20,4,22,10,22,12,20,12,16,10,14,6,12,6,8,50,50,5,4,7,4,7,6,5,6,5,4, 

/*  a  901. .942  */  9,18, 9,12, 8, 10, 4, 10, 2, 12, 2, 16, 4, 18, 8, 18,9, 16, 9, 12, 10, 10, 11, 10, 12, 12, 12, 18, 8,22,4,22,1 
/*  r  943. .950  V  8,22,4,22,4,4,8,4, 

/*  \  951. .954  */  1,22,11,4, 

/*  ]  955.-962  */  8,22,12,22,12,4,8,4, 

/*  A  963. .968  */  2,16,6,22,10,16, 

/*  _  969.. 972  */  0,4,10,4, 

/*  ”  973.. 976  */  4,22,8,16, 


/*  font  1001  symbols:  */ 


/*  000  */ 
/*  001  */ 
/*  002  */ 
/*  003  */ 
/*  004  */ 
/*  005  */ 
/*  006  */ 
/*  007  */ 
/*  010  */ 
/*  011  */ 
/*  012  */ 


3,2,15,20,3,20,15,2.3,2, 

3,11,13,11,50,50,8,7,7,8,8,9,9,8,8,7,50,50,8,13,7,14,8,15,9,14,8,13, 

3.6.15.6.50.50.15.8.3.13.15.18, 

3.6.15.6.50.50.3.8.15.13.3.18, 

4,3,14,3,9,14.4,3, 

2.14.4.9.7.19.17.19, 

6,8,12,14,50,50,6,14,12,8, 

3,17,15,17,50,50,3,11,15,11,50,50,3,5,15,5, 

2.14.14.14.50.50.2.8.14.8.50.50.4.4.12.19, 

2.3.3.2.5.2.6.3.10.19.11.20.13.20.14.19, 
4,5,12,5,14,6,15,8,16,11,15,14,14,16,12,17,4,17, 


/*  013  */  4,17,9,10,14,17, 

/*  014  */  4,10,9,17,14,10, 

/*  015  */  2,14,4,15,6,15,12,13.15,13,16,14,50,50,2,8,4,9,6,9,12,7,14,7.16,8, 
/•  016  */  14,4,14,3,3,3,7,11,3,18,14,18,14,17, 

/*  017  */  8,4,8,13,5,13,7,16,8,20,9,16,11,13,8,13, 

/*  020  •/  8,20,8,11,11,11,9,7,8,4,7,7,5,11,8,11, 

/*  021  */  2,11,10,11,10,14,14,12,17,11,14,10,10,8,10,11, 
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/*  022  */  8,11,8,8,4,10,1,11,4,12,8,14,8,11,16,11, 

I*  023  */  8,2,11,20,50,50,9,17,7,15,6,12,6,9,7,6,10,5,12,6,14,9,14,12,13,15,11,17,9,17, 

/*  024  */  14,12,5,12,4,8,4,6,5,4,7,3,9,3,11,4,12,6,13,8,14,12,14,16,13,18,12,19,10,19,9,18,7,16,5,12, 
/*  025  */  8,3,8,18,50,50,2,17,3,18,13,18,14,17, 

/*  026  */  8,2,11,20,50,50,3,16,4,16,5,15,6,12,7,10,9,9,11,9,12,10,13,12,14,15,14,16, 

/*  027  */  6,15,4,12,3,9,3,6,4,4,6,3,8,3,9,4,10,7,10,10,9,11,8,5,10,3,12,3,15,5,16,9,16,13,15,15, 

/*  030  */  3,19,5,19,6,18,12,4,14,3,15,3,50,50,9,11,5,4,3,3, 

/*  031  */  14,14,13,13,6,13,4,11,3,8,4,5,8,4,11,4,13,7,13,10,12,12,9,13, 

/*  032  */  13,19,12,20,10,19,8,17,8,14,9,11,11,10,12,7,11,4,9,3,6,3,5,5,5,8,7,10,9,11, 

/*  033  */  3,2,5,5,7,15,50,50,5,5,6,4,10,3,12,4,13,6,15,15,50,50,12,4,12,3, 

/*  034  */  4,3,6,12,50,50,11,3,13,12,50,50,4,11,5,12,7,12,9,11,10,11,13,12,15,12,16,13, 

/*  035  */  3,20,16,20, 

/*  036  */  4,8,14,8,50,50,4,15,14,15,50,50,9,10,9,20, 

/*  037  */  9,18,10,16,11,15,13,15,14,16,15,18,14,20,13,21,11,21,10,20,9,18, 


/*  font  0  lower  case  alphabet  characters  */ 


/* 

a 

1561. 

.1586 

*/ 

3,14,4,14,10,14,11,12,11,5,12,4,11,5,9,4,4, 

4. 2, 6, 2, 

10,3,11,11,11, 

/* 

b 

1587. 

.1606 

*/ 

2, 20,2,4,2,5,3,4,10,4,11,5,11,11,10,12,3,12 

,2,10, 

/* 

c 

1607. 

.1624 

*/ 

11,10,11,12,9,13,3,13,2,10,2,6,4,4,9,4,11,5 

, 

/* 

d 

1625. 

.1644 

*/ 

11,11,9,13,3,13,2,10,2,6,4,4,9,4,11,5,11,4, 

11,20, 

/* 

e 

1645. 

.1664 

*/ 

2,9,11,9,11,11,9,13,3,13,2,11,2,5,4,4,10,4, 

11,5, 

/* 

f 

1665. 

.1682 

*/ 

10,18,10,19,9,20,5,20,4,19,4,4,50,50,2,13,6 

,13, 

/* 

9 

1683. 

.1708 

V 

11,11,10,13,3,13,2,11,2,6,3,5,10,5,11,6,11, 

11,11,0, 

10, -1,3, -1,2,0 

/* 

h 

1709. 

.1722 

V 

2,20,2,4,2,12,3,13,10,13,11,11,11,4, 

/* 

i 

1723. 

.1744 

*/ 

6,4,8,4,7,4,7,12,6,12,50,50,6,16,6,17,7,17, 

7,16,6,16, 

/* 

j 

1745. 

.1768 

V 

4,1,5,0,6,0,7,1,7,12,6,12,50,50,6,16,6,17,7 

,17,7,16 

,6,16, 

/* 

k 

1769. 

.1780 

*/ 

2,20,2,4,2,9,7,15,2,9,10,4, 

/* 

l 

1781. 

.1790 

*/ 

6,20,7,20,7,4.6,4,8,4, 

/* 

m 

1791. 

.1814 

*/ 

1,13,1,4,1,11,2,13,6,13,7,12,7,4,7,12,8,13, 

10,13,11 

,12,11,4, 

/* 

n 

1815. 

.1828 

*/ 

3,13,3,4,3,12,4,13,10,13,11,12,11,4, 

/* 

0 

1829. 

.1846 

*/ 

11,6,11,11,9,13,3,13,2,11,2,6,4,4,10,4,11,6 

f 

/* 

P 

1847. 

.1866 

V 

3, -2, 3, 13, 3, 11, 4, 13, 10, 13, 11, 11, 11, 6, 10, 4, 4 

,4,3,6, 

/* 

q 

1867. 

.1886 

*/ 

11,-2,11,13,11,11,10,13,5,13,4,11,4,6,5,4,10,4,11,6 

1 

/* 

r 

1887. 

.1898 

*/ 

5,4,5,14,5,12,7,13,11,14,12,12, 

/* 

s 

1899. 

.1924 

*/ 

10,11,9,13,4,13,2,10,2,9,3,8,7,8,10,7,11,6, 

11,5,10, 

4, 4, 4, 2, 5, 

/* 

t 

1925. 

.1940 

*/ 

7,20,7,6,8,4,9,4,10,6,50,50,4,13,10,13, 

/* 

u 

1941. 

.1954 

*/ 

3,13,3,6,4,4,10,4,11,6,11,13,11,4, 

/* 

V 

1955. 

.1960 

*/ 

3,13,7,4,11,13, 

/* 

w 

1961. 

.1974 

*/ 

2,13,5,4,7,13,50,50,6,13,8,4,11,13, 

/* 

X 

1975. 

.1984 

V 

3,13,10,4,50,50,3,4,10,13, 

/* 

y 

1985. 

.2004 

*/ 

3,13,3,7,6,4,10,4,11,6,11,13,11,0,9,-1,6,-1 

,4,0, 

/* 

z 

2005. 

.2112 

*/ 

3,12,10,12,3,4,10,4 

>; 


void  gennote  (Twidth,  Theight,  Tlocate,  Trotate,  Tstring) 
float  Twidth,  Theight,  TlocateC3],  Trotate; 
char  ‘Tstring; 

C 

float  charwidth,  H_unit,  W_unit; 

Coord  xO,  xl,  yO,  yl ,  xt,  yt,  -C4]; 

register  int  i,  j; 

int  'Skip,  iptr,  nchars,  moveflg; 


nchars  =  strlen(Tstring); 
charwidth  =  Twidth  /  ( f loat Knchars) ; 
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H_unic  =  Theight  /  22.0; 

U  unit  =  charwidth  /  16.0; 


if  (  fabs(Trotate)  >  0.01  )  C 
/*  Text  is  to  be  printed  at  an  angle.  */ 


r [0]  =  cos(Trotate  *  M_P1 
r[1]  =  sin(Trotate  *  M_PI 
rC2]  =  TlocateCO]  *  (  1  - 
rC3]  =  TlocateCI]  *  {  1  - 


/  180.0); 

/  180.0); 

r CO)  )  +  TlocateCI)  *  r[1); 
r  CO)  )  -  T  locate  CO)  *  rCD; 


for  < i =0;  i<nchars;  i++)  < 
moveflg  =  0; 
iskip  =  0; 

iptr  =  *(Tstring++); 

for  ( j=(chstart Ciptr] - 1 );  j<=(chstopCiptr) -3);  j  +=  2)  C 
if  (iskip  ==  1)  { 
iskip  =  0; 
continue; 

> 

if  (charset Ci+2)  >  47)  C 
moveflg  ■  0; 
iskip  =  1; 
continue; 


> 

xO  =  TlocateCO]  +  (float)(  charset Cj]  )  *  W_unit; 
yO  =  TlocateCI)  +  (float)(  charsetCj+1)  )  *  H_unit; 
xl  =  TlocateCO]  +  (float)(  charsetCj+2]  )  *  U_unit; 
yl  s  TlocateCI)  ♦  (float)(  charsetCj+3]  )  *  H_unit; 


if  (  fabs(Trotate)  >  0-01  )  < 

/*  Apply  rotation/ translation  to  text.  */ 
xt  =  xO  *  r CO]  •  yO  *  rCD  +  rC2); 

yt  =  xO  *  rCD  +  yO  *  r  CO]  +  r  C3] ; 

xO  =  xt; 
yO  =  yt; 

xt  =  xl  *  r  CO]  •  yl  *  rdl  ♦  rC2); 

yt  =  xl  *  r  Cl)  +  yl  *  r  CO]  ♦  rC3]; 

xl  =  xt; 
yl  =  yt; 

> 

if  (Imoveflg)  move  (xO,  yO,  0.0); 
draw  (xl,  yl,  0.0); 
movef l g  =  1 ; 

> 

TlocateCO]  ♦=  charwidth; 

) 

> 


fortran  fcharst( i .TextStr,  n,  x,  y,  z,  Rotation) 
int  i; 

char  *TextStr; 
int  *n; 

float  *x,  *y,  *Z,  ‘Rotation; 
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extern  int  letterwidth,  letterheight,  strokewidth; 

float  Twidth,  Theight,  locationt3); 

char  s[80]; 

memcpy(s,TextStr,*n); 

s[*n]  =  *  \0  • ; 

locationCO]  =  *x; 

locationtl]  =  *y; 

location [2]  =  *z; 

Twidth  =  (<float)< letterwidth  '  strlenls))); 
Theight  =  ((floatXletterheight)); 
l inewidth(strokewidth); 

gennote< Twidth,  Theight,  location,  'Rotation,  s); 

linewidth(l); 

return(O); 


long  StrUidth(TextStr) 
char  'TextStr; 

C 

extern  int  letterwidth;  /*  letterwidth  is  set  in  Font(n)  */ 
return  ( long) ((strlen(TextStr))* letterwidth); 


fortran  f font( FontNumber,  BoldOpt) 
int  'FontNumber; 
char  8oldOpt; 

extern  int  letterwidth,  letterheight,  strokewidth; 
int  i;. 

i  =  'FontNumber; 
strokewidth  =  1; 
switch  (  i  )  < 

case  10:  letterwidth  a  10; 

letterheight  =  13; 
break; 

case  12:  letterwidth  =  12; 

letterheight  =  15; 

if  (  BoldOpt  ==  'b'  )  strokewidth  =  2; 
break; 

case  14:  letterwidth  =  15; 

letterheight  =  18; 

if  (  BoldOpt  ==  • b’  )  strokewidth  =  2; 
break; 

case  16:  letterwidth  =  16; 

letterheight  =  20; 

if  (  BoldOpt  ==  'b1  )  strokewidth  =  2; 
break; 

case  24;  letterwidth  =  25; 

letterheight  =  31; 

if  (  BoldOpt  ==  • b '  )  strokewidth  =  3; 
break; 

case  32:  letterwidth  =  30; 

letterheight  =  40; 

if  (  8oldOpt  ==  ■ b 1  )  strokewidth  =  3; 
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break; 

case  48:  tetterwidth  =  50; 

letterheight  =  63; 
if  (  BoldOpt  ==  1 b* 
break; 

case  8: 
default: 

letterwidth  =  7; 
lutterheight  =  10; 
break; 

> 

return(O); 

> 

long  GetHeightO 

return  ( long) letterheight; 

> 


)  strokewidth  =  4; 
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